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A Method For Numerical Integration 


By C. B. Haselgrove 


1. Introduction. In this paper we shall give an account of some methods de- 
veloped for the numerical evaluation of multidimensional integrals. These methods 
are based on the theory of Diophantine approximation. They are suitable for some 
problems for which the Monte Carlo method is commonly used and, like the Monte 
Carlo method, are well fitted for use with an electronic digital computer. We shall 
show, however, that they are superior to the Monte Carlo method provided that 
the integrand satisfies certain conditions. We shall also show that they are superior, 
for integrals in space of several dimensions, to formulas typified by those of Gauss 
and Simpson; they may be superior even to certain new integration formulas 
specially constructed for the evaluation of multiple integrals (see for example 
Hammer [2], who gives a bibliography, and Miller [5], [6], [7]). 

The method of antithetic variates which is described by Hammersley and others 
[3], [4] may be used to obtain better estimates than the Monte Carlo method but 
the author thinks that the method described in the present paper is simpler to apply 
and gives better results. 

Various authors have suggested methods which are particular cases of those 
described in this paper but without the underlying theory. See for example Davis 
and Rabinowitz [1]. 

In this section we shall give a short account of the behavior of the error in the 
Monte Carlo method and the direct-product Gauss-type methods so that we can 
compare these with the errors of the new methods. We shall not give an account of 
the method of antithetic variates. 

Suppose that we wish to estimate the integral 


1 1 1 
1=[ [ es [ flaiyaey +++ ye) dz dae +++ day. 
0 0 0 


We shall denote the vector (2, 2, --- , x) by x. Numerical methods for the 
evaluation of J involve the calculation of f(x) at a number N of points x; . The 
most desirable of such methods for use on an electronic computer are those which 
require the evaluation of f(x) at the smallest number of points x; in order to obtain 
an estimate with a sufficiently small error. 

The Monte Carlo method gives as an estimate for J the sum 


1 N 
7 2d S(xs) 


where the points x; are chosen at random in the range of integration. The error 
of such an estimate has standard deviation O(N~“’) provided that the function 
f(x) satisfies certain conditions. It is sufficient that the function be bounded. 

A Gauss-type formula for a one-dimensional integral takes the form 


a+h ” 
| g(y) dy =h > ¢;9(a + ah) + R. 
a i=] 
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The remainder R is 0(h’”) as h varies for functions g(y) satisfying certain condi- 
tions. We may estimate the integral 


ajith pagth agth 
/ / oe [ f(x) dz, diz see dx, 
@) ae aE 


by applying the formula in each dimension, obtaining the direct product estimate 


WD ++ Do ei +++ Cia + ah, +++, ae + ah). 
ij=1 ip=l 
This gives the value of the integral with an error which is 0(h‘**”) as h varies, 
provided that the function f(x) possesses partial derivatives of order up to 2» with 
respect to each variable and that these derivatives are bounded. These conditions 
are very strong in that f(x) must have the derivative 


Orkf 
onto? ---oa’ 
but they are probably more than is necessary. 
If we divide the hypercube 0 < 2; < 1 into n* smaller cubes with side h = 1/n, 
then using the »-point Gauss formula we obtain a formula for J which requires 
N = (m)* values of f(x). This gives the value of J with an error which is 


O(n-*’*") nin o(ns?”*?"*) 


for fixed v and k. Thus if 2v — 1 < $k, v and k being fixed, the error will not de- 
crease as N increases so rapidly as in the Monte Carlo method. This suggests that 
v should be chosen large. However, it is not desirable to use a high-order Gauss 
formula since large errors may then arise if high-order partial derivatives of the 
function are large or if the function does not possess such derivatives. The high- 
order partial derivatives are large if the function has singularities in the 2k-dimen- 
sional complex domain near the region of integration. 

Similar arguments to those above may be applied to multidimensional integra- 
tion formulas obtained as the direct product of rules such as Simpson’s rule, which 
gives an error that is O(N“). They may in fact be applied to any case where 
the domain of integration is divided into smaller regions as above. 

In this connection the following two integration formulas should be mentioned. 


@) app [2+ [1 aan --- dey = (1 - $100, 0, ---, 0) 


+ 4{f(h, 0, ---,0) + f(—h, 0, ---,0) + f(0, h, 0, ---,0) 


+ f(0, —h, 0, --+,0) + +++} + OCH). 


wie aml ‘@ [ 4) dey «+ dey = 37(0,0, ---, 0) 


h 


1 


+ aay [fh h, +++, h) + f(—h,h, ---,h) + flh, —h, h, ---,h) 


+ f(—h, —h, h, soo, h) + -+-} +0). 
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Formula (i) requires (2k + 1) values of f(x) to estimate the integral over the 
hypercube of side 2h, and formula (ii) requires (2* + 1) values. However, if the 
hypercube 0 S 2; < 1 is divided into n* smaller ones of side 2h then to evaluate 
I by using (i) requires (k + 1)n* + O(n**) values of f(x) and by using (ii) re- 
quires 2n* + O(n*~”) values. Simpson’s rule used in this way requires (k + 2)n* + 
O(n*) values. Thus, for fixed k the formulas (i) and (ii) and Simpson’s rule give 
errors which are O(N“), but (ii) is the most efficient. 

In this paper we shall describe methods which give errors which are O(N‘), 
O(N~*), O(N~*) and O(N~“*) when applied in any number of dimensions to fune- 
tions which satisfy appropriate conditions. The points x; at which the values of 
the function f(x) are required can be calculated by a computer more easily than 
the random numbers required for the Monte Carlo method. No coefficients, such 
as are required in the case of the Gauss-type methods, need be calculated, nor is it 
necessary to decide in advance the number of points to be used. 

In the following sections we first describe the methods as applied to the integra- 
tion of functions periodic in each of the variables. In Section 3 we describe calcula- 
tions leading to numerical estimates of some of the errors and to suitable sets of 
points x; at which to calculate the function. In Section 4 we describe the practical 
application of the methods to non-periodic functions, and in Section 5 we give an 
example. 


2. The Integration of Periodic Functions. 


2.1. General Theory. Let f(z: , x2, --- x1) = f(x) be a periodic function with 
period 27 in each of the k variables x; , z2 , --- , 2 . We shall describe a number of 
methods for estimating the integral 


tao | [ [sandr de. 


These estimates will be based on sums of the type 


(1) 8(N) = > enmf(2aemay , 2xmex , --- , 2xmax), 


where the a; are certain linearly independent irrational numbers and the cy,, are 
certain coefficients chosen so that s(N) — I as N — . We shall be particularly 
concerned with sets of coefficients cy», such that the number of non-zero coefficients 
in the sum (1) is finite and of order N. We shall give estimates for the difference 
s(N) — I in terms of N and of bounds for the derivatives of f(x). 

Certain formulas which are particularly simple to use are based on sums s( NV) 
which can be expressed in terms of the repeated sums 


Si(N) 


E se2eme 


and 


SN) = > Spa(m) 


m=( 
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where r 2 2. The S,(N) are modified Césaro means of the sequence f(21me), terms 
with positive and negative m being taken together. We shall show that the means 


1 T 
3(N) - oN + i Si(N) 
and 
—_ 1 T 
&2(N) bed (NV + 1): 1? So(N VY 


which are of the form (1), give good estimates of J for large N. However, although 
for functions satisfying certain conditions the estimate for the error s.(N) — J is 
better than the estimate for s,(N) — J, we find that the higher-order Césaro means 
do not give better estimates of J. We shall show that if f(x) satisfies certain condi- 
tions then the means 


l iQ T Y T 
= ——___________ { §.(2N = 
«(N) = ar pqpen pH (S2N + 1) — 280N)} 
and 
1 Y (OAT y T 
s,(N) = (N + 1) {S,(2N ) — 48,(N — 1)} 


give substantially improved estimates of the integral. Under these conditions we 
shall show that for large N, 


s,(N) — I = O(N”) 


forr = 1, 2,3 and 4. 


In the theory of the general method we are led to the consideration of the 
function 


ky(@) — > Cuat™, 


in terms of which we shall give estimates of the error. If for a particular set of 
coefficients Cym , 
K 
k < -——{_— 
| kw(6)| 3 | N sin 36 |" 

for some fixed constant K = K(r) and for all N and @, we shall say that the method 
has order r with constant K. For the sums s,(N), s:(N), s:(N) and s(N) the 
corresponding functions ky(@) are 


(1) ~ 1 sin (N + 3)0 
kw (0) = (N+1) sno ’ 





rp, _ 1 __ sin 4(N + 1)0 
O- W inte 


(9) = 1 sin’ 3(N + 1)@sin (N + 3)0 





~ (N + 1)2(2N + 3) sin? 46 
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and 





(4) ae 1 sin‘ 3(N + 1)0 
ks () (N+1)* sint}o - 
Hence these sums give methods of order 1, 2, 3 and 4 respectively. We shall prove 


that a method of order r gives an error which is O(N”) provided that the inte- 
grand satisfies certain conditions. 


A method may have order r for several different values of r. In particular, if 


—m?/2N2 


Cum = Coe , 


where the constant cy is chosen so that 
> Cum = :. 


it follows from the theory of elliptic modular functions that the method is of order 
r for all positive r. However, this method has certain disadvantages in practice. 


2.2. The Numbers a;. We shall now prove a lemma which will enable us to 
prove that there exist sets of irrational numbers a, which may be used to obtain 
good estimates of integrals. This lemma is a standard result in the theory of Dio- 
phantine approximation but we include the proof here for completeness. The lemma 
also shows that for most sets of numbers a; the error is not substantially worse 
than in the estimates given. We show in Section 3 how suitable sets of numbers 
a; may be found. 

We shall denote a set of integers (nm; , m2, --- , m) by n and we shall write the 
scalar product (mya; + mea, + --- + nox) in the form n-a. 

LemMa. Let @(n) be any positive function such that 


> 1/¢(n) = 1. 
Then there exist irrational numbers a, , a2, ++: , a, such that 
(2) (1 )o(me2) --- O(m,)|n-a — n| S 1/(k + 1) 


for all sets of integers m , m2, --- , Mm, n not all zero. Further, the measure of the set 
of a, a2,°°*, a withO FS a; S 1 such that 


(3) Ibd o(n1)o(n2) --- o(m)|m-a — n| S 6/(k + 1) 


is less than 6. 
We may suppose without loss of generality that 0 < a; < 1. We observe that 
for fixed n, n the measure of the set of @ such that 


(4) (b(n) --- o(m)| n-@ — n| S 6/(k + 1) 
is less than or equal to 


6 1 1 





(k + 1) max | ni | o(m1)@(M2) -- + O(n) © 


But n is a non-negative integer less than or equal to k max | n; | . Hence the measure 
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of the set of a which satisfy (4) for given n and some n is less than or equal to 


5/(m1)o(n2) --- o(nr). 
Thus the measure of the set of @ which satisfy (3) is less than 
io} 2° 1 
3 fey = 8 
x x $(1)b(M2) «++ o(nx) 


(we have strict inequality since the regions overlap). Taking 6 = 1 we obtain the 
first part of the lemma. 





2.3. Estimates for the Error in the Integration of Periodic Functions. We now 
consider the problem of estimating the integral of a periodic function of k variables 
which can be represented by an absolutely convergent multiple Fourier series 

f(x) = +++ ¥ anjya,..-n,e™*. 
m1 Nk 
We shall suppose that for some positive fixed number s there exists a vonstant 
M, such that if none of the n; = 0, 


8 


(5) | Gnyng---ny | SS M,| myn --- m| ™. 


We suppose further that if any of the n; = 0 the same inequality holds with the 
zero factors omitted from the denominator on the right hand side. 

Since the Fourier series is absolutely convergent we may treat the contributions 
to the sum s(N) from the individual terms of the Fourier series separately. We 
obtain 


s(N) = >> ewmf(2ema) 


Il 


D +++ Do dnyng---nykw(2an- a). 
m1 


Nk 


Now 4o,.,--... = J, which is the integral which we wish to evaluate. We shall 
suppose that the coefficients cy, are scaled so that ky(0) = 1. Then 


(6) 8(N) = I = DS dnjng---nkw(2an- aa), 


where the prime ’ means that the term (0, 0, --- , 0) is omitted. 
We suppose that the method is of order r so that there exists a constant K such 
that 


for all N and 6. Now 





To 1 or 
N sin 36 | 


| sin ré | > 2 || € || 


where || é || denotes the distance of — from the nearest integer. Thus if the numbers 
a; satisfy the condition (2) we have 


|sin mn-a| = 2 l ’ 
k + 1 (mm) (m2) --- (nx) 


Thus if we apply the bound for the Fourier coefficient and the condition (7) we 
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deduce that the contribution of the term m , m,--- , m to the sum is jess than 
KM 2" (k + 1)"{o(m)o(me) --- o(me)}" | mime --- me | “N~. 


We have supposed that we may drop any of the n; from the denominator if it 
is zero. If we suppose that the sum 


> ‘to(n)}" | n | 


converges we deduce that there exists a constant C,,, depending only on r and s 
such that 


|s(N) —I| S$ KM.27(k + 1)'C,N. 
We have the bound 


Cre S {9(0)}’ + Di'fe(n)}" | n|~ 


for the coefficients C,,,. If s > r + 1 we can choose an appropriate function ¢(n) 
which leads to a numerical bound for the error of the method. However, this bound 
for C,,, is by no means the best possible and we shall show how it may be improved. 

We shall now prove the theorem. 

THEoreM. If the numbers a, a2,--- , ag are chosen appropriately and if the 
function f(a , 22, + , Xe) ts periodic with period 2x in each of the variables x; and 
its Fourier coefficients satisfy the condition (5) then if we apply any method of order 
r < 8, where r = 1, with constant K there exist numbers C,., such that 


(8) | s(N) — IT | S KM.27(k + 1)"ArstCeN~. 


The numbers C,,, depend only on the numbers a; and not on the particular method or 
the function to be integrated. 

If r = s it can be proved that | s(N) — J| = O(N***) for « > 0, but we 
shall not give the proof here. 

In order to prove the theorem we suppose that the function ¢(m) is chosen to 
be monotonic increasing for n = 0 and such that ¢(—n) = ¢(n). We divide the 
range of summation of the variables m , m2, --- , m into zones defined by 


N; < Nn; <= 2N; 


(or if n; is negative — 2N; < n; < — N; ; zones for which some of the variables are 
zero are also allowed). Now if n’ and n” are two sets of numbers in the same zone 
we deduce that |n,’ — n,” | = N;. Hence 


IV 


\|(a’ — n”)-@ || 2 1/(k + 1)6(m’ — m”)p(m2’ — m2”) --- (me’ — me”) 
1/(k + 1)6(N1)6(Ne2) -- + 6(Ne) 


IV 


so that 
| | n’-@ || — || n”-a|| | = 1/(k + 1)(Ni)o(N2) --- o(Ne)- 


We deduce that the number of sets of numbers n in the zone with 


|| m-a@ || S »/(k + 1)@(Ni)o( Ne) --- oC Ne) 
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is at most 2v + 2. Now for any set of numbers n, 
|| m-a@ || = 1/(k + 1)6(m)b(me) --- o(nx) 
so that the contribution of that set to the sum (6) is less than or equal to 
27° (k + 1)"M.K{o(m)ob(ne) --- b(me)}" | mune --- | “N. 


If we assume that {¢(n)}" || ~“* decreases monotonically for n > 0 we deduce 
that the contribution of this term to the sum is less than 


2° (k + 1)"M.K{o(Ni)o(N2) --- o(Nx)}" | NiN2--- Ni | “N. 


Thus, the sum of the contributions of all terms in the zone N; < n; S 2N; to the 
sum (6) is less than 


2°(k + 1)'M, Kto(Ni)o(N2) --- o(Nx)} Palla 
-|NiN2 +++ N,|7* (2 a ae ") N~. 


y=l 


Similar results hold for those zones for which some of the n; are negative or zero. 


If s > r and r > 1 we may choose ¢(7) so that the sum over all zones con- 
verges. We obtain the inequality 


|s(N) —1| 
< 2(k + 1)"M.K{2¢(r) + 2}({o(0)}" + 220 fo(ND} NT)“, 

where N, runs through the powers of 2 and ¢(r) is the Riemann zeta function. 
This inequality is of the form (8) with A, = {2¢(r) + 2}. 

If r = 1 and s > 1 we obtain an expression with a different C,,, and a coefficient 
A,,»,x depending on s and k. 

These estimates for the errors are very crude but give some indication of the 
power of the methods. In the next section we shall describe how sets of numbers 


a; ean be obtained and how numerical bounds for the errors in the estimates of the 
integrals can be calculated. 


3. Numerical Estimates for the Error. 


3.1. A Method for Obtaining Numerical Estimates. In this section we shall be 
especially concerned with the methods using the sums s.(N) and s(N) to estimate 
the integral of a periodic function 

f(x) — f(a eee "Ss 9 r) = > a Zz. @ayng---agt 
ny Nk 
Note that this function has period 2 in each of the variables xz; . As before, we shall 
assume that the Fourier coefficients satisfy an inequality 


(9) } Baishins + ing | s M, | nynz e's m| 


We shall be particularly concerned with the cases s = 2 and s = 4. In order to 
obtain a convenient form for the final results we shall replace zero factors in the 
denominator of the right-hand side of (9) by (6/2") when s = 2 and by (360/7z") 
when s = 4. 
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We recall that the functions ky(@) for r = 2 and r = 4 take the form 


(2) = 1 sin’ 3(N + 1)@ 
kw (@) = ap aint Fo 





and 


: 1 sin‘ 4(N + 1)@ 
kw (6) = (N+1)* sint3o ~ 





Both these functions are positive. 
Thus 


|s(N) —T| = DS) --+ Do “aajng---ngk” (am-@) 
M.>, --->.| mm --- m| ks? (an-a). 
We now define the functions 


f.(x) = M,* >, eee »™ | 24 Ne ++) Mh pinta 


k 20 
_ M,* II ( > infer) 


j=l n=—o 


lA 


where M,,* is chosen so that f,(0,0, --- ,0) = 1 and where, when n = 0, the fac- 
tor n * is replaced by the appropriate fraction for s = 2 or s = 4. 
By summation of the trigonometric series we deduce that 


P 
fo(x) = Tia —|2;|)9} 


and 
k 
f(x) = []{2— G — Jas) — | a )* 


We denote the values of the sums s,(N) corresponding to the functions f,(x) 
by s,,.(), and the mean value of f,(x) over the periodic cell by J, . 
Then we deduce that 


| s-(N) ae I | ~ (M,/M,*){8,,.(N) Tr T,): 


Thus, we may say that the functions f.(x) and f,(x) are the worst functions for 
numerical integration whose coefficients satisfy the inequality (9) with s = 2 and 
s = 4, 

If, therefore, we find sets of numbers a; which make (s,,,(N) — J,) for a prac- 
tical range of N as small as possible, these a; will give good estimates for the in- 
tegral of any periodic function whose Fourier coefficients satisfy (9). 


3.2. The Determination of the a;. The Ferranti Mercury computer at Man- 
chester University was used to find good sets of a; for (r,s) = (2, 2) and (2, 4). 
This was done by minimizing the upper bounds of (N + 1)*{s,,(N) — J,} for 
0 < N S N,. The minimization was performed by a random walk in the e-space, 
a step being taken only if it decreased the upper bound. The size of the step was 
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reduced when the upper bound could not be decreased otherwise and the number 
N, was increased when the step was reduced. 
For s = 2,7, = (4)* and for s = 4,7, = (4£)*. In both cases 


f(0,0,---,0) =1 


and f(x) 2 0 for all x, so that S,(m) = 1 and S.(N) = N + 1. Hence, for s = 2 
we have 


(N + 1)*{S22(N) — Io} 2 (N +1) — (N + 1)7(4)*. 


For a certain value N,, of N the right-hand side of this inequality approaches its 
maximum value of }-3*. Therefore, for values of N,; > N,, we cannot expect the 
upper bound lub (N + 1)*{8:2(N) — J:} to be less than 3-3". The corresponding 
value for s = 4 is 3(45)*. 

It was found that the points a converged to points for which the upper bound 
was reasonably small. In Figure 1 the logarithms of the upper bounds for s = 2 
and s = 4 are plotted against k. For s = 2, N, was taken up to 1500; for s = 4, 
up to 1000. The values }-3* for s = 2 and 3(4)* for s = 4 are shown on the same 
logarithmic graph as straight lines. The point for s = 2, k = 8 lies below the ap- 
propriate line; for this case N; was less than the number N,, . The numbers a; ob- 
tained for s = 2 and s = 4 are given in Tables 1 and 2 respectively. It should be 
realized that although these are good sets of a; they are not claimed to be the best 
possible. 


The convergence of the points a suggests that there may be limit points e’ such 
that the bounds 


lub (N + 1)"fs,,(N) — I} 
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Fig. 1.—Errors associated with the numbers a; for (r, s) = (2, 2) and (2, 4). The points 
are the logarithms to base 10 of the computed errors and the lines are those below which, for 
large enough N, the points cannot be expected to fall. 
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TABLE 1 


METHOD FOR NUMERICAL INTEGRATION 


Numbers a; for (r, 8) = (2, 2) 











0.73258893 0.95734608 0.80638723 
we 0.86730270 0. 22584927 
0.09724025 0.72510075 
0.62055505 0.31301950 0.51310685 
0. 22610245 0.48476582 0. 11080509 
=3 aoe 0.60161858 
0.92715171 
0.96498949 0.43657951 
0.81091316 0.59185199 | k=8 
0.46960090 0.05024400 0.73750248 
bin 0.84373919 0.08314415 
0.38104000 | 0.84753682 
0.62366851 | 0.75808683 | 0.88989711 
0.04150108 
| 0.80254484 
0.48574769 0.27951501 
0.27210703 0.67340402 
| 0. 53040927 





for these points are not substantially larger than the values givén in the figure, but 
the author has not been able to prove this. 

In the case k = 1,r = 2, s = 4, however, he has been able to prove the following 
result. 

If, for some a, (N + 1)*{se4(N) — I} forO < N < N, then there exists an 
irrational a’ such that for all N 


(N + 1)*{84(N) — Ij S 2Bi1 + O(N7*)} 
and such that | a’ — a| = O(Nz'™*). 


4. Application to the Integration of Non-Periodic Functions. We shall now show 
how the methods of Section 2 may be applied to obtain estimates of integrals of 
non-periodic functions of several variables. Suppose that it is required to evaluate 
the integral 


I= [/ cof Plz, a0, +++ 24) drs da --+ de 
R 


where R is some region of the k-dimensional space and F(x) is some given function. 
We shall show how under certain special circumstances it is possible to transform 
the problem to that of the integration of a periodic function. 

We have seen that we can obtain good estimates for the error in the integration 
of a periodic function if its Fourier coefficients tend to zero rapidly. Now, the 
Fourier coefficients tend to zero rapidly if the function is continuous and possesses 
partial derivatives of high orders. We shall, therefore, attempt to construct 
“smooth” periodic functions. The methods given below for doing this are intended 
only as examples. It will be apparent that each problem may require special treat- 
ment to obtain the best results. Also, if the function F(x) possesses unknown dis- 
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TABLE 2 
Numbers a; for (r,s) = (2, 4) 

0.83969144 0. 44810200 0.58505729 
k=2 0.53589831 0.50196855 
r 0 .56039410 | 0.77797734 
0.59734470 0.83630131 0. 60504620 
0.92828094 0. 22148205 0.62193588 
k=3 k=6 | 0.84244165 
“aahaaare 
0.74235492 0. 10613747 | simmers 

0.57387033 0.40278232 | k= 8 
0.32279917 0.88772556 | 0.23975940 
oars 0. 43554826 | 0.01544979 
0.17219381 0.57794809 
0. een 0.63794472 0.81182909 
= 3 tine 0.78068912 
-98875 | 0.62319488 
0.60299793 0.70710061 








0.60389317 


continuities, it will not be possible to remove these from the periodic function 
constructed. 
We first consider the integral 


_ 1 
I = | [  [ F(a, 22, 72+) tp) dx, day -+- dx,. 
0 0 0 


This integral may be written in the form 


1 1 1 1 : 
r=] [++ [ Pet lel +5 |eel) den dee --> dep. 
—ji oe] i 


The function f(z , 22, --: , 2) = F(\a\|,|2\|,---,|2|) may be regarded as a 
periodic function with period 2 in each of the variables z;. Then f is continuous if 
F is continuous, and if F possesses a derivative 


gt 
O21 0X2 -** OXy 


of bounded variation it may be proved that the Fourier coefficients of f satisfy, 
for some M; , 


| Qnyng---mg | S&S Me | nyne --- me |. 


This implies that the means s,(N) and s:(N) will give estimates for J with errors 
O(N") and O(N~***) respectively. These methods have the particular advantage 
that since f(z) = f(—z) the means s,(N) and s:(N) may be evaluated from the 
values of F at N + 1 points. The sum S,(N) may be expressed in the form 


N 
Si\(N) = F(0,0,--- ,0) + 22 F(2 \{3nay}| , -- + , 2 |{3naz} |) 
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where the brackets { | denote the fractional part, lying in the range (—43, 3). For 
such a function the basic points of Table 1 are suitable. 
We now consider the application of higher-order methods to the integral 


1 esl 1 
I = | il tee [ F(z, 22, +++, 2) dey dzz --+ dz. 
—} 1 1 
We make the substitution z; = P,(z;) where r = 1 and 
P,(2;) = A, [ (1 — wv) du 
0 
and A, is chosen so that P,(+1) = +1. Then 


I 


1 1 
[ me? [ Fea, +++, P(a,))P,/(x4) sa P,!( ax) dx, --- dx, 


Il 


1 1 1 1 
I [ ves [flay xa, +++, 24) day dry +-- da, 
= af 1 on 


say. Then f = 0 if any z; = +1, and F is bounded. Thus if F is continuous f may 


be regarded as a continuous periodic function of the z;. If F possesses a partial 
derivative 


Eoeer Mas: 
(d2, O22 ay = Ox;)" 
of bounded variation then so will f. It is then possible to obtain an estimate 
s(N) —I = O(N”) 


so that the application of a method of order r will give an error of order N” and 
the application of a method of order r + 1 will give an error of order | 

The numbers a; of Table 2 can be used to calculate the sums s.(N) when the 
partial derivative a” f/(Ax0x2 --- d2,)* of the function has bounded variation. 

It is also possible to apply these methods to the more difficult problem of the 
evaluation of an integral taken over some more general regions. It is, of course, 
possible that there may exist some transformation of the variables which reduces 
the integral to periodic form. For example, in. the evaluation of the integral of a 
function taken over the interior of a circle we may take polar coordinates. The 
function is then a periodic function of the angular variable. It may be made into 
a periodic function of the radius with derivatives of appropriate orders by one of 
the techniques above. 

In the case of the integration of a function F(z) over the infinite range (— ~, ~ ) 
we may at once consider the periodic function 


= +e 





> F(z + 2xn) 


—2 


but this may mean that we have to take a very large number of points if F(z) does 
not tend to zero rapidly ast > +. 

A more general method is the method of extrusion. Suppose that a region R with 
a smooth boundary is contained in a sphere 8 with center O in R. Then, if P is 
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any point in the sphere we may construct the line OP. If R is a star body about O 
this line will meet the boundary in one point Q. We may then evaluate the function 
F at a point P’ at a distance OP-OQ/r from O along the line OP, where r is the 
radius of the sphere. If this value is multiplied by (OQ/r)* we obtain a function 
defined in the sphere which has the same integral over the sphere as F has over R. 

This process of extrusion may be generalized to the extrusion of any star body 
to any other star body containing it. It will generally be better to use a sphere 
rather than a cube if R has a smooth boundary, otherwise discontinuities are in- 
troduced into the derivatives. 


5. Example of Use. The method with the most practical application is probably 
that using the sum s,(N). The sum s,(N) may be more suitable if higher accuracy 
is required but for the results to be better than those for s.(N), the function f must 
have higher-order partial derivatives. 

The values of s,(N) can be calculated at chosen values of N while the repeated 
sum S,(N) is being accumulated. The accuracy of the results can be estimated by 
examining the convergence of the successive values of s:(N). 

A refinement may be introduced in order to decrease the rounding errors. The 
sum S2(N) becomes very large and for a floating point computer with a fixed 
number of significant figures the rounding error introduced at each addition is 
proportional to the larger of the two quantities added. Thus, the sum s.(N) will 
not be accurate to the full capacity of the computer. A method of overcoming this 
difficulty is to calculate an estimate J,* = S.(N)/(N + 1)’ at intervals in the 
computation, and to subtract this value from each subsequent calculation of the 
integrand, so that the accumulated sums are 

N 


Si*(N) = Dd {f(me) — 1,4 














m=—N 
and 
N 
S.*(N) = > Si*(m). 
TABLE 3 
Estimaies for fi fi eee fi errryr"ns dx, dx, eee dxs 
N s2(N) | a(N) Monte Carlo 

1000 0.97062580 |  0.97062392 | 0.96763166 

2000 0.97063927 0.97082902 0.96870265 

3000 0.97066765 0.97054070 0.96885258 

4000 0.97066383 0.97068153 0.96944396 

5000 0.97065630 0.97065925 0.96950137 

6000 0.97065761 0.97061983 0.96990269 

7000 0.97065639 0.97068925 0.97018578 

8000 0.97065632 | 0.97064881 0.97030504 

9000 0.97065706 0.97063833 0.97038771 
10000 0.97065854 0.97066307 0.97032729 
11000 0.97065860 0.97065947 6 .97029480 
12000 0.97065744 0.97067426 Cc 








Exact value = 0.97065719 
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The sum S,*(N) does not then become very large. At each successive re-estimate 
Tixs = 1,* + Ss*(N)/(N + 1)? we replace the current S,*(N) by zero and the 
current S,*(N) by S,*(N) + (2N + 1)(1,* — I%4;:). The same method used in a 
fixed-point computer will reduce the amount of scaling necessary. 

As an example the integral 


[f ee [ ett de da, -- 


was calculated on the Mercury computer. Some results for k = 5 are given in Table 
3. We give the sum s.(N) for N = 1000( 1000) 12000. The method described above 
was used to reduce rounding errors, with the estimate J,* calculated at every 10th 
value of N up to N = 100 and every 100th value thereafter. For comparison, the 
sum 3,(N) is also given, and an estimate using the ordinary Monte Carlo method 
with N random sampling points. 
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The Fairing of Ship Lines on a High-Speed 
Computer 


By Feodor Theilheimer and William Starkweather 


Abstract. Methods for using a digital high-speed computer to determine ship lines 
are presented. It is assumed that the offsets of a small number of points were taken 
from a preliminary design, and that it is desired to compute the offsets of an arbi- 
trarily large number of points on the ship’s surface. Procedures for using a computer 
for the solution of this problem are described. Special emphasis is placed on the 
detection, by a computational criterion, of unwanted fluctuations and the correc- 
tion of such fluctuations if they should occur. The method also includes a special 
procedure which takes care that those portions which are straight in the preliminary 


design remain straight in the final form. Illustrative examples of the methods are 
discussed. 


1. Introduction. Before the actual construction of a ship can begin, considerable 
time and effort has to be spent in the drawing of the ship lines. The purely graphical 
methods of determining ship lines are very tedious and time-consuming. Therefore, 
the problem of implementing the graphical methods by analytic procedures has 
been studied for a considerable time. 

One of the most important earlier contributions is due to Admiral David W. 
Taylor [1]. An extensive history and bibliography of the problem is given in Volume 
II. of “Hydrodynamics in Ship Design,” by Captain Harold E. Saunders [2]. 
Among the more recent publications are papers by W. H. Résingh and J. Berghius 
[3], P. C. Pien [4], and J. E. Kerwin [5]. 

The analytical approach to the ship line problem has become particularly at- 
tractive since high-speed computers became available. This offers an opportunity 
to eliminate much of the drudgery inherent in the graphical method. 

When ship lines are found by an analytic method, they may possess unwanted 
fluctuations. It is, therefore; desirable to have an analytic criterion which permits 
us to determine whether or not a line is free of unwanted fluctuations. This paper 
furnishes such a criterion and gives a method of finding lines without such fluctua- 
tions. 

In Section 2 a method, which essentially amounts to an interpolation, is de- 
veloped for finding a ship surface which passes exactly through a set of points 
given in a preliminary design. 

In Section 3 this method is modified to a smoothing procedure which yields 
ship lines that are free of unwanted fluctuations but which may no longer pass 
exactly through the given points. 

Section 4 deals with the special case where the preliminary design contains 
straight line portions, and a method which reproduces these lines as straight lines 
is developed. 


Received May 15, 1961. This paper is reprinted from David Taylor Model Basin Report 
1474, January 1961. 
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7 
Fic. 1.—The Coordinate Axes 


Section 5 describes applications of the various procedures for finding ship lines 
by a high-speed computer. 


2. Interpolation Method. The coordinate system is taken as indicated in Figure 1; 
the z-axis is longitudinal, the z-axis is vertical, and the y-axis is perpendicular to 
both the z- and z-axes. 

The determination of ship lines by analytic methods is considered equivalent to 
finding a function of two variables 


(1) y = F(z, z) 


so that for every fixed value of z we have y as a function of z along a waterline, 
which permits us to compute an arbitrarily large number of half-breadths on that 
waterline. Likewise, a fixed value of x would yield an arbitrarily large number of 
half-breadths along a section. 

We assume that the ship for which we are to determine y = F(z, z) is described 
to us by a table of offsets taken from a preliminary design. 

Let the offsets be given at N + 1 points on each of K + 1 waterlines. Then 
we will develop a procedure which yields y = F(z, z) and thus gives an arbitrarily 
large number of offsets instead of the given (N + 1)(K + 1) points. The existing 
program for the computing machine is set up for handling cases where N and K 
can be as large as 24, but, in the cases actually treated, N and, particularly, K 
were much smaller. 

To find a function y = F(z, z) we first determine a function y = f (x) which 
corresponds to a single waterline. The finding of such a function y = f(z) will 
turn out to be an essential part of the determination of the surface y = F(z, z). 

We assume that the entire waterline is divided into N, not necessarily equal, 
intervals and that the N + 1 points (2, yo), --- , (tw, yw) are given. We now 
consider the problem of finding a curve which passes exactly through these points. 
This interpolation problem can be attacked only after we make a choice as to the 
type of curve that should represent a waterline. 

We shall insist that the function be continuous and have continuous slope and 
curvature or, what amounts to the same, have continuous first and second deriva- 
tives. 

A simple type of curve that satisfies these conditions is the following: 

Let the curve consist of a number of segments, each segment being represented 
by a cubic. These cubics are joined in such a way that at the juncture points the 
function, its first derivative, and its second derivative are continuous. 

Curves of this typé arise in the theory of small deflections of thin beams which 
are simply supported at a finite number of points. To some extent a batten or 
spline held in place by so-called ducks, as it is used in the drawing of ship lines, 
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can be approximated by a thin beam simply supported at a finite number of points. 
The analogy between a spline and a thin beam gave rise to the name “spline curve” 
for the type of curve consisting of cubic segments described here. 

We have chosen spline curves to represent ship lines. The explicit formula for 
a spline curve is greatly simplified if we introduce the following notation: 


lA 
~ 


0 for x 


(z — a),’ 


= (x — a)’ for x 


IV 
& 


Equation (2) can also be written in the form: 


(3) (x — a)y* = 3[(x — a)’ +|z—a| 
With the aid of this notation, a spline curve which has discontinuities of the 
third derivative at the points 2, , --- , zy, can be written as 


f(z) =y=at+br+ cr + Ayr’ 
+ A(x = nm) + --- + Ayu(« = 2u-s)+ 


It is clear that each new term 


(4) 


A,(x — n)4°, n=1,---,N-1 


introduces a discontinuity of the third derivative at x = x, , the magnitude of the 
jump in the third derivative being equal to 6A, , whereas the continuity of the 
function, its first derivative, and its second derivative remain undisturbed. 

In each of the N intervals 


Gena < 2X Bus n=1,2,---,N 
formula (4) sums to a single cubic. In individual cases any of the four terms of such 
a cubic may drop out in the summation, and, in particular, it is possible that the 
cubic may reduce to a straight line. This shows that the spline curve lends itself 
to the representation of curves which contain straight portions as, for instance, 
waterlines on ships with parallel middlebody. 

We shall now develop a procedure for determining the coefficients of a spline 
curve which passes exactly through the N + 1 points 


(2o , Yo), (x1 ) yi); “a cil (ty ’ yw). 


We choose as points of discontinuity of the third derivative the N — 1 inner 
points 


(2 ’ ¥); ine (ty ) Yn-1) 


of the given set of points. The spline curve then appears in the form of equation (4) 
which has N + 3 coefficients. For the determination of these coefficients we have 
N + 1 equations; namely, the conditions that the curve should pass through the 
N + 1 given points, y, = f(z.) forn = 0,1, --- , N, which means that we have 
two more unknowns than equations. 

To get a clearer picture of the role of these two degrees of freedom, we rewrite 
f(x) of equation (4) in the form 


(5) f(x) = fo(x) + pfi(x) + afe(zx) 
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where fy , f: , and fz are to satisfy the following conditions: 


-o( Zn) = Siltn) _ 0 fo(2n) _ 0, eS = 0, 1, ° > N 
(6) fo'(a) = 0 fi' (ao) = 1 fe! (x) = 0 
fo" (x) = 0 fi’ (xz) = 0 So” (x) = 2. 


Each one of the functions fy , f; , fe is now determined by N + 3 conditions. The 
requirements at 2» are clearly satisfied if we write: 


fo(x) = Yo + B(x > I)” 
+ Biz — %)4° + --- + Bya(z — tr-4)+? 


(7) f(z) = 2 —2y + C(x — %)* 
é 
+ Ci(x — )+ + --- + Cyi(r — Zs) 4° 
f(z) = (x - 2)” + Dox — 2)” 


+ D(x — m)4° + --- + Duala — ty) 4°. 


The B, , C, , and D, ,n = 0, --- , N — 1, can be found by utilizing the equa- 
tions fo(tn) = Yn, filan) = 0, fo(an) = 0,n = 1, --- , N. To find these coefficients 
one has three systems each of N equations in N unknowns, which are obtained by 
substituting the values of x; , 22, --- , 2» in place of z in the three equations of 
(7), and by utilizing the first line of equation (6). The coefficients B, , C, . and 
D, can be found without solving the systems of equations in the usual manner. 
Substituting x = 2, immediately yields B, , Cy , and Dy , and if we put z = z, , we 
find B, in terms of B,, B,,--- , B,-, and similarly, C, and D, in terms of the 
preceding C’s and D’s, respectively. 

After the functions f, , f: , and fe are thus determined, the parameters p and ¢ 
of equation (5) must be found. Among the curves 


y = fot ph + gfe 


which all pass through the N + 1 points (2, yo), --- , (tw , yw), we have to choose 
one which is by some criterion most desirable. 

Between two spline curves, we will consider as more desirable the one which 
has smaller jumps of the third derivative, or more precisely, the one where the 
sum of the squares of the jumps of the third derivative is less. Going back to thin 
beam theory, we note that the reaction force at a point of support is proportional 
to the jump in the third derivative. The desired spline curve is, therefore, the one 
for which 


(8) » (Ax)? 
becomes a minimum, since 6A, is the jump in the third derivative at x, for the 


function f(x) of equation (4). 
From equations (5) and (7) we see that 


(9) A, = B, + pC,r + QD, 
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and therefore 


N-1 
(10) >» (B, + pC, + qD,)° 
n=l 
has to become a minimum. 
This leads to the following two linear equations in p and q 
p2c,, + qzC,D,, - — B,C, 
p=C,.D, + q=D,” = — B,D, 


(11) 


where we always sum in n from 1 to N — 1. 
By this method a spline curve which passes through N + 1 given points is 
found in the form: 


N-1 
(12) f(x) = yo + p(x — x) + g(a — 20)? + Ag(x — a)* + >. A,(z — 2n)a°. 
n=l 
By expanding the powers of (x — 2»), this can be brought into the form 
9 n—1 " 
(4) f(x) = a + br + cx + Agr’ + D> An(t — 2n),°. 
n=l 


We have thus developed a procedure for finding a spline curve which represents a 
waterline that passes through N + 1 given points. If we are given N + 1 points 
on K + 1 waterlines, we can therefore find K + 1 functions, f,(2), k = 0,--- , K, 
expressing the K + 1 waterlines. Our next step will then be to find y = F(z, z) 
so that 


(13) F(a, 2) = f(x) 


where 2; is the z coordinate of the kth waterline. 

Each one of the functions f,(2) appears in the form of f(x) as given in equation 
(4) where the coefficients a, b,c, Aj, Ai, --- , Aw-1 are generally different for each 
f(x). Each one of these coefficients, say for example, a is now given for z = 2, 
2=%,°°',2 = 2c. We can therefore ask for a function a(z) for which we are 
given values at 2, 21, °°: , 2x. This is exactly the same problem that confronted 
us when we were given N + 1 points on a waterline and found f(z). 

By using the same procedure which gave us f(x) we can now find a(z) as a spline 
function of z. In the same manner b(z), c(z), Ao(z), --- , Aw-1(z) can be found. 
This gives the expression for the whole surface 


F(a, 2) = a(z) + b(z)-a + c(z)-x” + A(z) -2° 
(14) 


+ A,(z)-(a — 1)4° + ---+ Ay-i(z)-(2 = Sys). 


Each one of the N + 3 coefficient functions a(z), b(z), e(z), Ao(z), A1(z), 
Ay-i(z) is itself a sum of K + 3 terms. 

In this derivation we worked under the assumption that the abscissas 2, , 
n= 1,---, N —1 do not change from waterline to waterline. In general, we 
assume that the data given originally will satisfy this requirement. However, we 
permit 2) and zy , which do not appear explicitly in the formula for the spline curve, 
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“ 


to vary from waterline to waterline. This is very advantageous since, in general, 
waterlines on different levels do not have the same z, and zy . 

In order to apply our procedure as it stands, it is necessary that z,,--- , 2w-4 
do not change from waterline to waterline. However, a simple modification of our 
procedure will permit us to handle the case where some waterlines do not conform 
to this restriction. 

If some waterlines have offsets given at 2, 2, --- , Z» , whereas one waterline 
has offsets given at 2*, 2:*, --- , y*, then we can find a function which expresses 
this waterline in a form analogous to equation (13) by means of a function 


f(z) = a+ br + cx” + Agr’ + As(z — a*)4° + °°: + Aya(z — 2y-1)4° 


From this function we can find f(x), f(ae), --- , f(@w-1) which can serve as ordi- 
nates of modified data which go with the abscissas x, , --- , Zy-; . Such a modifi- 
cation can also be used if one line contains fewer points than the other lines. 

The derivation of F(x, z) described here shows immediately that the originally 
given points on each waterline do not have to have equidistant abscissas nor do the 
waterlines have to be equally spaced. 

A program to compute the coefficients of the function y = F(z, z) was coded for 
the Remington Rand high-speed computer UNIVAC. Further computer programs 
were then developed to evaluate this function and to print out in tabular form the 
offsets along any portion of any waterline if z is given and the first desired value of 
x, an increment of x, and the number of z values desired are entered into the ma- 
chine. In similar fashion, any portion of any transverse section can be tabulated. 

These various machine programs were later recoded for the IBM 704 computing 
system. 

It is clear that the determination of the coefficients for the function y = F(z, z) 
constitutes the main computing effort. The evaluation of the function and the 
printing of the results for any number of points requires very little extra time. 


3. Smoothing Method. The procedure described thus far usually gives very satis- 
factory results, but in some cases it turned out that the computed ship lines showed 
unwanted fluctuations which could be recognized by plotting the lines. There 
arises a need for a criterion which will permit the computer to determine whether 
or not a line is free of such unwanted fluctuations. 

Clearly, a curve that has no inflection points would be free of all fluctuations. 
We cannot require, however, that ship lines have no inflection points, because 
there are well-designed ship lines that do have inflection points. What we want to 
avoid are extraneous inflection points. We therefore have to determine where the 
given data indicate the presence or absence of an inflection point. 

This can be done by studying three successive points P,.. , P,» , Pas: with the 
coordinates (2n—1, Yn—-i), (Xn, Yn), (Xn4i, Ynsi)- If the three points are plotted and 
P,,_, and P,, 4; are joined by a straight line, it can readily be seen whether the curve 
in the vicinity of P, is concave upward or downward. This graphical method can be 
easily translated into analytical form by introducing the second difference, which 
here shall be called r,, . 

Then the sign of the second difference will indicate whether the curve is to be 
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concave upward or downward at P,,. A positive second difference indicates con- 
cave upward, a negative second difference indicates concave downward. 
The formula for the second difference becomes particularly simple in the case 
of equally spaced abscissas: 
(2n41 —_ Z_) = (Za wa Zn—1) = h. 


Then the second difference becomes 





fe Yn et 2Yn + Yn-1 
(15) t, = 73 ‘ 


In the case where the abscissas are not equally spaced, formula (15) is to be re- 
placed by the general formula 


(16) ee 2 Ynti — Yn _ Yn — Bo) 
n Tat — Ln-1 \Cnti — Ta 4. — Ss > 


If we are given the N + 1 points 





(Zo, Yo), ee (tw , Yw) 


we can readily determine N — 1 second differences 
"1, “ee » nw-1- 


If r, and r,4; have the same sign, we can consider this an indication that there 
should not be an inflection point between z, ard 2,4: . If they have different signs, 
the given data suggest an inflection point between z, and z,,;. We are now in a 
position to determine whether or not a curve y = f(z) has unwanted fluctuations. 
Only if the second derivative f”(z,) and the second difference r, have the same 
sign can we be sure that the curve y = f(z) is free of unwanted fluctuations. 

Here it is understood that a function whose second derivatives at zx, and 2n4; 
have the same sign has no inflection points between z, and z,4; . This, however, 
will not be true for every function; for instance, it will not be true for a polynomial 
of degree higher than three. Such a function can have an even number of inflection 
points between two successive points at which the second derivatives have equal 
signs. However, if we deal with spline curves, z, and z,4; are joined by a cubic, 
which can possess not more than one single inflection point. By the same reasoning, 
we find that if in the case of a spline curve, f”(z,) and f”(x,4:) have different signs, 
there will be precisely one inflection point in the interval z, to 2,4; , and not an 
arbitrary odd number as it might be for the general curve. 

For illustra**on, we have discussed the comparison of the second differences of 
the data with the second derivatives of the computed curve only in the case where 
the curve has been a waterline, and we have, therefore, been finding a function of zx. 
Exactly the same test can, of course, be carried out when we are finding a function 
of z. 

The test of whether or not a given curve has unwanted fluctuations can be 
carried out very easily by the computing machine for either waterlines or sections, 
and the program can be written so that the computing machine prepares a listing 
of the points where the curve fails the test. It may occasionally happen that the 
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discrepancies are of a trivial nature, but if they are not, a method must be devised 
to compute new lines which are free of unwanted fluctuations. 

In trying to free ourselves from unwanted fluctuations we think of a smoothing 
procedure which relaxes the requirement that the computed curve should pass 
exactly through the given points. 

The method of least squares immediately suggests itself. If, however, a least- 
squares method is applied, one finds that in comparatively simple cases one is led 
to unwanted fluctuations. The fact that a least-squares fit may yield fluctuations 
can be ascribed to the following consideration. 

A least-squares fit minimizes the sum of the squares of the differences between 
the given and computed values but imposes no constraint to minimize fluctuations. 
If we want to free ourselves of fluctuations, we have to take care that the second 
difference and the second derivative have the same sign in each case. This can be 
attempted by minimizing the sum of the squares of the discrepancies between the 
second differences and the second derivatives. 

Such a minimization is essentially implicit in the interpolation procedure de- 
scribed by equations (4) through (10). There we minimize (A, )*. It can be shown 
that 


(17) An = 5 lta — $"(a0)] 


in the case where the abscissas are equally spaced with the spacing A. 
Let z, be a point where the third derivative of a spline curve undergoes a jump 


(18) 6An = f"(an4) — f"(t0-) 

then 

(19) flan +h) = flats) + hf'(ae) + *$" (20) + <I" (aus) 
(20) fle — h) = flan) — hft(an) + % g(a.) — % mx, 


If we add equations (19) and (20) and utilize equations (15) and (18), we 
get equation (17). 

Therefore, when we used the interpolation’ method described by equations (4) 
through (10), we found among all possible spline curves passing exactly through 
the N + 1 given points the one for which, in the case of equally spaced points, the 
discrepancies between the second differences and the second derivatives were least. 

Even if we minimize the discrepancies between the second differences and the 
second derivatives but insist that the curve pass exactly through the given points, 
we will often get unwanted fluctuations. To free ourselves from these fluctuations, 
we have to relax the requirement that the computed curve pass exactly through the 
given points. 

If we no longer insist that the curve pass exactly through the given points, we 
do not require that f(z,) be equal to y, , but we still want to make 


u [f(xn) 9 Yul 
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small. To avoid extraneous fluctuations we also want to make >—*=) [f’(2,.) — ral 


n=1 


small. We therefore minimize a linear combination of the two sums; namely, 
N N-1 

(21) Do [f(an) — yal + 820 If" (an) — rol 
n=0 n=l 


The use of the parameter s enables us to emphasize either the first or the second 
sum. By choosing s small we put more stress on minimizing the first sum, which 
makes the curve pass more closely to the given points. By choosing s large we put 
more stress on minimizing the second sum, which makes the discrepancy between 
the second derivative and the second difference small. This tends to make them 
agree in sign, and thus tends to eliminate unwanted fluctuations. 

If we study the case of s = 0 we minimize only the first sum. This sum has a 
zero minimum, since each summand can be made zero individually by solving the 
N + 1 linearly independent equations, f(z.) = y., n = 0, 1,---, N, for the 
N + 3 unknown coefficients, a, b, c, Ag, -'- , An-1. Therefore s = 0 leads to a 
curve passing exactly through the given points. This curve, however, would not be 
uniquely determined but would depend on two parameters just as in equation (5). 

A similar situation would prevail if we were to minimize only the second sum in 
formula (21), which would be equivalent to letting s become infinite. Then we would 
also get a zero minimum, and we would achieve perfect agreement between the 
second derivatives and the second differences. For sufficiently large s we can there- 
fore make the second sum as small as we please and obtain agreement in sign be- 
tween second derivatives and second differences. 

Thus we have found a method for computing a curve which is suggested by 
N + 1 points. We first minimize expression (21) with a small s and then test to 
see if there are any unwanted fluctuations. If there are some, we increase s until 
the curve is satisfactory. When all waterlines have been fitted with satisfactory 
curves, we can then fit like coefficients into functions of z by utilizing formula (21) 
again, and thereby find the expression for the ship’s surface y = F(z, z). 

Although one may choose different s values for different lines and even for dif- 
ferent points on a given line, most of the work has been done primarily by using 
the same value of s throughout for any given run on the computer. 


4. Lines With Straight Portions. One of the reasons for choosing spline curves as 
typical ship lines was the fact that ship lines often contain straight portions. A 
polynomial or any other analytic function cannot contain straight as well as curved 
portions. Although spline curves can contain straight portions as well as curved 
ones, the methods developed thus far do not contain any special provisions which 
will ensure that straight portions will be reproduced as straight portions. 

It is felt that every effort should be made to reproduce exactly those portions 
of the preliminary design which are clearly intended to be straight. 

First, we will treat the case where the straight portion is at the beginning of 
the line. 

Suppose that of the N + 1 points 


(Zn » Yn), n = 0,::- »N 


the first L + 1 points lie on a straight line of slope a. 
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Then 


(22) Yn = Yo t+ A(X, — 2X) forn = 0,---,L 
and the desired function can be written 
N-1 
(23) f(x) = yo + a(x — x) + 2d Aalz — 2) 4°. 
The N — L unknown coefficients A, ,n = L, --- , N — 1 could be found by 
solving the N — L equations: 
f(2r41) = Yrsr = Yo t+ a(2r41 — to) + Ar(tins — 21)’ 
S(tr42) = Yrs2 = Yo t+ a(2142 — Zo) 
+ A1(4142 - tz)° + Arai(Ti42 = Zi43)° 
(24) 
Yw = Yo + a(tw — X%) + Ar(tw — 21x) 
+ Ary(tw — Zi) +--+ + Aya(tw — tw)’. 


After having done this we would have to apply the previously derived criterion 
to find out whether or not there were any fluctuations. If there were unwanted fluc- 
tuations, we would then minimize the expression 


f (xn) 


N N-1 
(25) D (fan) — yo +8 DO [Ff (an) — rol. 
n=L+1 n=L+1 

As before, a sufficiently large s will ensure that the sign of the second derivative 
f” (xn) of the computed curve will agree with the sign of the second difference r, of 
the given data at each z, . 

If a ship line consists of a curved part followed by a straight part of slope b; 
that is, if, of the N + 1 points, 


(Zn, Yn); n=0,---,N 
the last L + 1 points lie on a straight line, then 
(26) Yn = yw — O(ty — 2.) 
forn = N,N —1,---,N-—L. 


The expression for the ship line is then 


1 
(27) f(z) = yw — b(ayw — x) + De An*(tn — z),’ 
where (z, — x),° = Ofor (xz, — x) < 0. 

The A,,* can be determined by a procedure analogous to the one used in equa- 
tions (24) or (25). The expression for a line thus determined will contain terms of 
the form A,*(xz, — xz),°, whereas we usually write lines with terms of the form 
A,(x — &na)4°. 

The two forms are related to each other in the following way 


(28) (x — an)4° — (an — 2) 4° = (2 — 2p)’. 
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That this is a valid relationship becomes evident when we realize that for z less 
than z, , the first term on the left-hand side of equation (28) becomes zero, whereas 
the second term vanishes for x greater than z, . 

Therefore, equation (27) can be modified in the following way: 


1 


f(x) —— b(ty — z)+ ae A,*(a, — z),’ 
(29) n=N— 


1 1 
= E —Way—z)-— DY An*(2 - r.)'| + DO An*(2— 2n)4°. 
n=N—L n=N—L 
The terms inside the brackets can be arranged as a single cubic. 

We have thus far treated the cases where a straight portion is followed by a 
curved portion and where a curved portion is followed by a straight portion. There 
is one more possibility; namely, that a line has a straight portion followed by a 
curved portion followed by a straight portion. If there are one or possibly more 
straight portions in the interior of a line, then we can subdivide each of the inner 
straight portions and treat the thus achieved cases individually. 

The case of a curved portion between two straight portions leads to the problem 
of finding a spline curve for which the ordinates, the slopes, and the second deriva- 
tives of the two end points are given. In particular, the second derivatives are 
zero. Since six conditions are imposed, the spline curve must have at least six 
coefficients, which means that it must consist of at least three cubic arcs. 

We first consider the case where there are two points given between the end 
points of the straight portions. Let x, be the abscissa of the end point of one of the 
straight lines and z,4: be the initial abscissa of the next straight line. 

If 


f(t1) = A f(%i43) = B 
(30) f(t.) =a f'(z14s) = b 
f’ (a1) = 0 f” (2143) = 0 
then 
L+3 
(31) f(z) =A+ a(x — 2.) + D> Aa(z — 2,),° 
n=L 


which immediately satisfies the conditions at x, . 

The conditions at x,4; yield three linear equations in the three unknowns A, , 
Ars, and A,42. The determinant of this system is recognized to be different from 
zero because it contains three rows which are the first, second, and third powers of 
(ti4g — An), N= LE+1,L4+2. 

The coefficient A,4; is determined by the requirement that for x > 2,4; the 
curve should be a straight line; namely, f(x) = B + b(x — 2,43). 

The solution of this problem can be given explicitly for the case of equidistant 
abscissas ©, , X, + h, x, + 2h, x, + 3h. To simplify the necessary formulas, we 
look at it as a combination of three elementary cases. 

In the first case, if f(z.) = A and f(z, + 3h) = B and f'(z.) = f’(2z1) = 
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f'(zi + 3h) = f’ (2. + 3h) = O, then by solving equations (30) adapted for this 
case, 


(32) f(z) =A+ a (x — 21)+° — X2@-—2,- h),? 


+ Xz -2- 2h),? -—(t#-a- 3h) 4’). 





In the second case, if 


f(ar) = f(a + 3h) = f’ (x1) = f'(zi + 3h) = f" (zi + 3h) = 0, 
but f’(z.) = a, then 


(33) S(#) = a(z — 21) + So [-2( — 21)4? + 5(@ — 2, — hi), 
— 4(x — 2, — 2h)4’ + (x — x, — 3h), "I. 


In the third case, if f(z.) = f(r. + 3h) = f’(z1) = f’ (21) = f" (ai + 3h) = O, 
but f’(z. + 3h) = b, then 


(34) f(x) - aal-(2 om tr)4° + 4(a —-F4-— h),° 
— 5(z — 2, — 2h),* + 2(x — x, — 3h),'). 


If f(tz) = A, f(t. + 3h) = B, f'(z.) =a, f'(t. + 3h) = 6, f"(21) = 
f’ (x. + 3h) = 0, then f(z) can be written as the sum of the three expressions 
(32), (33), and (34). 

It is to be observed that for both the equidistant and the nonequidistant solu- 
tions of equations (30), the curved part is completely defined by the slopes and 
the end points of the given straight portions. The magnitudes of the ordinates at 
141 and 2,42 do not enter into the computation; therefore, if fewer than two 
points are given in the gap between the end points of the two straight portions this 
method can also be applied. All that is necessary is to choose additional abscissas 
in the gap without fixing their ordinates. 

A different situation prevails if the curved part between the two straight portions 
has more than two given points. Let the first straight line end at z, so that 


(35) f(t) =A, fi(%)=a, f(t.) =90 
and the next straight portion begin at zy so that 
(36) f(zu) = B, f'(zm) = b, f" (tm) = 


where M — L > 3. 

We can immediately give one such spline curve which satisfies at least all con- 
tinuity requirements at its junctures with the two straight portions. This can be 
done by extending the first straight line until only three intervals have to be bridged 
by the curve. The following relations will have to be satisfied: 


fo(am—s) = A + a(ru_s — Zz) fo(zu) = B 
(37) fo' (tus) = @ fu' (am) 
fo (rus) = 0 fo (rm) 
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0. 
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The conditions at z+ y_; are immediately satisfied when we write 


fox) = A + a(x — 21) + Aw_s(t — tu_s)4° + Auo(t — tu-2),° 
(38) 
+ Aw-a(z — Zu-i)+ + Ax(z — tu)+. 


The coefficients A y_3;, Am—2, Am_i follow from the conditions at xy whereas A y 
is a result of the requirement that the line should be straight for x > ry. 

Although such a spline curve satisfies the continuity requirements at x, and 
Xm, in general, it cannot be expected to pass closely enough to the interior points 
(Xr41, Yrs), (Ve42, Yr42), °** » (Wu-1, Yui), nor is there any assurance that the 
curve will be free of fluctuations. 

The particular solution of the problem of finding a spline curve for which equa- 
tions (35) and (36) hold can be expanded into a general solution if we add linear 
combinations of spline curves which satisfy the conditions 


(39) f(ai) = f'(ar) = f’ (aL) = f(au) = f'(am) = f’ (am) = 0. 


The general solution satisfying (39) is of the form 
M 
(40) f(z) = > A,(2 — 2)+. 


The first three conditions of equations (39) are automatically satisfied. The last 
three conditions of equations (39) lead to three equations in A,, Ar4y,---, 
Awu-i, Whereas Ay is determined from the requirement that the spline curve be 
straight for z > ru. 

The three equations in A,, --- , A w_; have coefficients which form a matrix of 
rank three, since it contains three rows which are, respectively, the third, second, 
and first powers of 

(tm —2,) forn=L,L+1,---,M —1. 


Since the difference between the number of unknowns and the rank of a linear 
equation system gives the number of linearly independent solutions, we will have 
here M — L — 3 linearly independent solutions. 

One way of finding a special set of such independent solutions is to determine 
M — L -- 3 functions which satisfy the conditions 


(41) fi(x) = 0 fort, St S tips 
and for tr4i43 S 2 S 4m 
and which in the interval r,4.-1 < © < 21443 is controlled by the conditions 


Se’ (@r44-1) = 0 Sa’ (Xi 4e43) = 0 
(42) 
fi (tr40-) = 0 Sr (214043) = 0 


fork = 1,2,---,M—L-—83. 

These equations (41) and (42) are a case of equations (39) if we replace L 
by L+k—1 and M by L + k + 3, and thus, M — L becomes 4. Therefore, 
for every value of k = 1, 2,---, M — L — 3, we get exactly one independent 
solution, fi.(x). 
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Since along with f,(z) any multiple of f,(z) is a solution of these equations, we 
normalize by requiring that 


(43) filZise) = 1. 


These functions f,(z) are then shown to be linearly independent in the following 
way. 
If there were to exist a relation 


M—L-3 


> Vihlr) =0 
k=l 


then substituting x = 2,4, would yield V, = 0, and, successively, all V; could be 
shown to be zero. 

The general solution of equations (35) and (36) is then found by adding to (37) 
the general solution of (39), or what amounts to the same, a linear combination 
of solutions of (41) and (42). 

We thus obtain 


(44) fix) = fox) + > pefe(2). 


The coefficients p, are to be determined in such a way that 


(45) ¥S Ue) —yk ts © (en) — rl 
n=L+1 n=L+1 


becomes a minimum. 

The larger the value of s the more emphasis is placed on the agreement of the 
second derivative f”(z,) of the curve with the second differences r, of the given 
data. It is to be observed that in this case we have M — L — 3 coefficients, and 
we want to have agreement between the second derivatives and the second dif- 
ferences at M — L — 1 points, so that no matter how large an s is taken we can 
only minimize the discrepancies but cannot force the discrepancies to zero. In 
most cases we will easily achieve agreement of sign between the second derivatives 
and the second differences with a suitable value of s. 

Should a case arise where it cannot be achieved, we can take recourse to the 
following procedure. Agreement in sign between the second derivative and second 
difference can be forced if we permit some of the quantities to vary that have been 
held constant previously. Such quantities are, for instance, the constants which 
characterize the straight portions; namely, the quantities a, A, b, B, of equations 
(35) and (36). 

It may be of interest to note that in the case of equidistant abscissas zx, , the 
solution of equations (41) and (42) can be given explicitly in comparatively 
simple form. 

If x14. = x, + kh, then a simple verification shows that 


g(x) = (a — 2,)* — 4(a — a, — h)4* + 6(2 — 2, — 2h),° 
(46) 
—_ A(x = = 3h), + (Z—- @- 4h).° 





352 FEODOR THEILHEIMER AND WILLIAM STARKWEATHER 


can be taken as f(x), and all the f,(z) are found by the formula 
(47) fi(x) = glx — (k — 1)hl. 


Thus, we have developed a procedure for finding spline curves into which we in- 
corporate straight portions. The methods described apply to abscissas that are not 
necessarily equally spaced, but they simplify considerably in the case of equi- 
distant abscissas. 


5. Applications. The various procedures described so far were applied to both 
realistic and artificially created data. We now want to describe in some detail two 
particular problems which were solved on the IBM 704. 

In the first case, a set of faired lines for DD 948-949 were available. A small 
number of offsets were selected as the basis of a problem. We applied our method 
to them as if they were taken as part of a preliminary design, and reconstructed 
from them a large number of points in order to compare them with the correspond- 
ing faired data. 

In particular, eighteen points on each of five waterlines were chosen. The chosen 
waterlines were 2 ft, 4 ft, 6 ft, 10 ft, and 14} ft. The point (zo , yo) on each water- 
line was the tabulated initial point of that line so that z» was not the same for all 
waterlines. The other abscissas 2, , --- , X17 were chosen so that z, = 20 n ft. 

With these 90 given points, the coefficients of the function representing the ma- 
jor part of the ship’s underwater surface were determined by the method of Sec- 
tion 3. 

This permitted us to compute any number of offsets and, in particular, we were 
able to tabulate the offsets at 2-ft spacing on each of the input waterlines and also 
on the 8-ft and 12-ft waterlines. 

These computed offsets were then compared with the given faired offsets and 
the discrepancies were found to be no greater than $-in. at any point. 

This example, used as a practice case, was based on input points which were 
already well faired; therefore, a small value of s, namely, s = 1, in formula (21) 
gave no fluctuations. Also, when the interpolation method described in Section 2 
was applied to the same data no fluctuations appeared. 

A practical application of the methods described here was made in the solution 
of a problem proposed by the New York Naval Shipyard. 

It was requested that faired ship lines for LPD 1 be computed. A preliminary 
design was given, including a table which contained the offsets on each of 11 water- 
lines for stations that were generally 25 ft apart. The problem was to compute the 
offsets for each frame at 2-ft spacing. 

A study of the preliminary design of the ship showed a feature which required 
some extension of the method described previously. Thus far we have considered 
only ships whose lines have continuous slope and curvature at every point. This 
ship, the LPD 1, has a horizontal knuckle 22 ft above the keel extending from the 
midsection to the after perpendicular. At this knuckle the section lines have a dis- 
continuity of the first derivative. 

The waterlines can be handled here in the same way as previously described. 
If we were to continue as in the case where there is no knuckle and form a surface 
y = F(a, z), we would be led either to violent fluctuations or, if s in formula (21) 
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TABLE 1 
Waterline Half-Breadths teat New York Naval Shipyard for 
LPD 1 
(From Station 18, z = 450 ft, to Aft Perpendicular, z = 500 ft; 
Station Spacing 25 ft. Dimensions in ft-in.-8ths) 














Station | 22-Fi W.L. | 28-Ft WL. 36-Fi W.L. | 4-Fi W.L. 
18 | 35- 1-4 | 35-11-4 37-0-6 | 38-2-0 
19 | 3111-4 32-11-6 344.2 | 35-90 
AP | (28 40 | 29- 7-0 313-4 | 33-00 





were chosen sufficiently large, the knuckle would be completely faired out, which 
is not what was desired. 

We can avoid this difficulty by subdividing the coefficients of the various water- 
lines into two groups: those that go with the fore part of the ship where there is no 
knuckle, and those that go with the aft part of the ship where there is a knuckle. 
The fore part is treated in the usual manner. In the aft part, we fit the coefficients 
of the waterlines above the knuckle separately and also the coefficients of the water- 
lines below the knuckle separately. 

In the actual computation, formula (21) was used with s =~ 1, which gave satis- 
factory results. The result of the computation was printed by the high-speed 
printer associated with the computing machine. It was arranged in the form of 
tables giving offsets on every frame. The frames were spaced 2 ft apart, and on 
each frame the offsets were given for every waterline level 2 ft apart. 

As an illustration, a small part of the input data furnished by the New York 
Naval Shipyard is presented in Table 1, and the corresponding output computed 
by the IBM 704 is given in Table 2. 

Table 1 contains the portions above the knuckle of Stations 18, 19, and the After 
Perpendicular. 

The computed data for the same region are given in Table 2 for Frames 225 
to After Perpendicular (Frame 250), which corresponds to a length of 50 ft, and for 
the 22-ft to 44-ft waterlines. 

Tables 1 and 2 contain what is essentially 1/20 of the actual input and output 
data. ' 


6. Conclusion. The methods derived here were developed to serve as a means for 
finding an arbitrarily large number of full-scale offsets if only a limited number of 
offsets, taken from a preliminary design, are given. This particular use guided us 
in the choice of the type of curve to represent a typical ship line. As was made 
clear, particularly in Section 3, special care was taken to avoid unwanted fluctua- 
tions. 

If we fit a whole waterline or a half waterline with a single polynomial, we do 
not have any assurance that our lines are free from unwanted fluctuations. Even 
a least-squares procedure will not necessarily exclude unwanted fluctuations, since 
it controls only the discrepancies between given and computed points and does 
not put any constraint on the curvature. In Section 3 we described a method where 
such a constraint on the curvature is applied. 
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TABLE 2 


Waterline Half-Breadths Computed on IBM 704 for LPD 1 


Dimensions in ft-in.-8ths) 


(From Frame 225, z = 450 ft, to Frame 250, z = 500 ft; Frame Spacing 2 ft. 





22-Ft W.L. 


























24-Ft W.L. 26-Ft W.L. 28-Ft W.L. | 30-Ft W.L. 32-Ft WL. 
225 | 35-14 35-4-7— | 35-81+ | 35-114 | 36-26+ | 36-61 
226 | 34-10-5+*| 35-21- | 35-54— | 35-87— | 36-02- | 36-3-5- 
a7 | 34-7-64 | 34-11-2 35-2-6— | 35-6-1+ | 35-9-5— | 36- 10+ 
22g | 34- 4-7 34-8-3+ | 3411-7+ | 35-33+ | 35-67+ | 35-10-34 
299 | 34-20-*| 34-544 | 34-9-1— | 35-054 | 35-4-2— | 35- 7-64 
230 | 33-110- | 34-25 34- 6-2 34-9-7 | 35- 1-4 35- 5-1 
231 | 33-80- | 33-11-5+ | 34-33- | 34-7-0 | 34-106- | 35-23 
a32 | 33-4-74+ | 33-85+ | 34-034 | 34-41 | 34-7-7 | 34-115 
233 | 33-1-7— | 33-55 33-9-3+ | 34-12 | 34- 50+ | 34- 8-7 
234 | 32-10-6— | 33-244 | 33-634 | 33-10-2+ | 34-214 | 34- 60+ 
935 | 32-7-4+ | 32114— | 33-33 33-7-3— | 33-11-2 34- 3-2— 
236 | 32-43-— | 32-83— | 33-0-2+ | 33-4.2+ | 33-82+ | 34 03- 
237 | 32- 1-1 32-5-1+ | 32-92— | 33-12 33- 5-3— | 33- 93+ 
238 | 31- 9-7 32-2.0-- | 32-61— | 3210-2 | 33-23- | 33-64 
239 | 31-6-5— | 31-106— | 32-27+ | 32-7-1— | 321124+ | 33- 34+ 
240 | 31- 3-2 3i-7-4— | 31-116— | 32-40- | 32-82 | 33-044 
241 | 30-11-7 31-4-1+ | 31-84— | 32-0-6+ | 32-5-14+ | 32-94 
42 | 30-84— | 31-0-7— | 31-5-2— | 31-95 32-204 | 32-6-4- 
243 | 30-5-0+ | 30-94— | 31-1-7+ | 31-63 3110-7 | 32- 33+ 
244 | 30-1-5— | 30-6-1— | 30-105— | 31-3-1 31-7-6- | 32-03- 
245 | 29-10-1— | 30-25 30-7-2— | 30-11-7 31-44 | 31-9-2- 
246 | 20-6-44+ | 29-11-14 | 30-3-7— | 30-84+ | 31-1-2+ | 31- 60+ 
247 | 29-3.0— | 29-7-5+ | 30-0-3+ | 30-5-2— | 30-10-0+ | 31-27 
248 | 28-11-3- | 29-41 29-9-0— | 30-1-7— | 30-66 | 30-11-54 
249 «| 28-7-54+ | 20-0-4+ | 29-5-4— | 29-10-34 | 30-33+ | 30- 8-4— 
250 | 28- 4.0 28-8-74+ | 20-1-7+ | 2-7-0 | 30-0-1— | 30- 5.2- 
Frame 34-Fi W.L. 36-Ft W.L 38-Ft W.L. 40-Fi W.L | 42-Fi W.L. | 44-Fi W.L. 
225 | 36-9-3+ | 37-06 37-4.0+ | 37-7-3- | 37-10-5+ | 38-20 
226 | 36-7-0— | 36-10-3- | 37-16— | 37-5-1— | 37-84— | 37-11-7- 
207 | 36-44— | 36-7-7+ | 36-113— | 37-26+ | 37-62 | 37- 95+ 
22g | 36-1-7+ | 36-5-3+ | 36-87+ | 37-04— | 37-40- | 37-74 
229 | 35-113— | 36-2-7+ | 36- 6-4 36-10-0+ | 37-1-5+ | 37-52 
230 | 35- 8-6 36- 0-3 36- 4-0 36-7-54+ | 36-112+ | 37-30- 
231 | 35-6-1— | 35-9-6+ | 36-14 36- 5-2— | 36-8-7+ | 37-0-5+ 
232 | 35-3-3+ | 35-7-1+ | 35-10-7+ | 36-26- | 36-64 | 36-103- 
233 | 35-0-5+ | 35- 4-4 35- 8-3 36-0-2— | 36-4-1— | 36- 8.0- 
234 | 34-9-74+ | 35-1-7— | 35- 5-6 35-9-5+ | 36-15 | 36- 5-5- 
235 | 34-7-14+ | 34-11-1 35- 3-1 35-7-1— | 35-11-1 36- 3-1+ 
236 | 34-43- | 34-83 35-0-3+ | 35-44 | 35-85- | 36- 0-6— 
237 | 34- 1-4 34- 5-5 34- 9-6 35-17 | 35-6-0+ | 35-10-2— 
238 | 33-10-5+ | 34-27- | 34- 70+ | 34112 | 35-3-4- | 35-7-6- 
239 - 76 34- 0-0 34- 4-24 | 34-85- | 35- 0-7 35- 5-2— 
240 | 33-4-7— | 33-9-14+ | 34-14 34-5-7 | 34-10-2 | 35- 25+ 
241 | 33-1-7+ | 33-6-2+ | 33-106 34-3-1+ | 34-7-5 | 35-0-1— 
m2 | 32-11-0— | 33-33+ | 33-7-7+ | 34-034 | 34-5.0- | 34- 9-4 
243 | 32-80— | 33- 0-4 33-5-1— | 33-9-54+ | 34-224 | 34- 67+ 
244 | 32-4-7+ | 32-9-5- | 33-29- | 33-67 33-11-5— | 34- 4-2 
245 | 32- 1-7 32- 6-5 32-11.3— | 33-4-1- | 33-87— | 34-145 
246 | 31-10-7— | 32- 3-5 32-8-4— | 33- 1-2 33- 6-1— | 33-10-7+ 
247 | 31-76 32- 0-5 32- 5-4 32-10-3+ | 33-3-3- | 33- 82— 
248 | 31- 4-5 31-9-5— | 32-25- | 32-74+ | 33-04+ | 33- 5-4 
249 | 31-14 31- 6-44+ | 31-11-5 32- 45+ | 32-9-6— | 33-246 
250 | 30-10-3- | 31-34— | 31-85 32- 1-6 32-6-7 | 

















* Plus sign indicates +, in.; minus sign indicates —,; in. 


33- 0-0-— 
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The representation of a single line by means of different cubics for different 
parts of the line, which might be cumbersome for hand calculations, does not con- 
stitute any difficulty for a high-speed computer. 
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Symmetric Successive Overrelaxation In Solving 
Diffusion Difference Equations 


By G. J. Habetler and E. L. Wachspress 


1. Introduction. In [1] Sheldon presented an iteration scheme for solving certain 
elliptic difference equations. The computational experiments described in his paper 
indicated that this method was superior to the method of successive overrelaxation 
[2]. We shall show that this conclusion is valid for the model problems that he con- 
sidered but that it is not valid for general diffusion difference equations. 

In this paper, theoretical convergence of the method is established, a variationai 
scheme for estimating the optimum parameters is developed and numerical results 
are given for some typical diffusion calculations encountered in nuclear reactor 
theory. 


2. General Theory. The two-dimensional diffusion difference equations can be 
simplified to be of the form 
(1) (I — B)o = S, 


where B isareal symmetric n X n matrix with spectral radius less than unity and 
with zero diagonal] entries, S is a known source vector, and ¢ is an unknown flux 
vector. The Jacobi method of iteration, which is also known as the method of simul- 
taneous displacements [2], for solving system (1) is given by 


(2) o” sini Bo“? + a. 
Let E“ equal @ — ¢'”, where ¢ is the unique solution to equation (1). The error 
vector E™ satisfies 
(3) E© = BE“ 

Consider some ordering ¢ of the integers 1 < i S n. Let P, stand for the cor- 
responding n X n permutation matrix and let 

B. = P.BP,’. 

Then B, can be written as B, = R, + R,” where R, is a lower triangular matrix 


with zero diagonal entries. Now let ¢, = P.@ and S, = P,S. The method of suc- 
cessive overrelaxation is given by 


(4) af" _ (1 = ow + wlRos” + DN + S,]. 
The error vector satisfies 
EY” = LewEs’ = (I — wR) {(1 — w)I + oR."JEs 


or in terms of the ordering implied in equation (1), 
E® aii P a F P EC 
(5) b . 
= (I — oR,)"[(1 — w)I + oR,"|E° 
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where &, = P,"R,P, . If the spectral radius of B (the maximum of the-magnitude 
of the eigenvalues of B) is 7(B), Young [2] has shown that for “consistent order- 
ings” the optimum value of w to use is 


2 
© T+ V1 — UB) 
For this value of w, the spectral radius of L,.,, is 
(6) a2.) 2 bs ME. 
: 1+ V1 + FB) 


If a(B) = 1 — &e<1), then alL,..] = 1 — 20/2e + (ce). If we define the 
rate of convergence of the successive overrelaxation method, R(L,...) as 








(7) R(Le.o) = —In aL.) 
we see that, for 7(B) close to unity, 
(8) R(Lg,0) & 20/2. 


The method of symmetric successive overrelaxation with extrapolation is a three- 
step process given by 


¢ < (1 ay wo? + wR?” 4p Rg? + S] 
9) aa an (1 ad wae + lke”? + R76 + S] 
or = is + af¢?” ay ” in + ia” a g*) 


where the a, and b, are chosen according to a technique developed by Stiefel [3], 
an extension of a polynomial extrapolation scheme introduced by Shortley and 
Weller [4]. The error vector for this scheme satisfies 


(10) E® = MS2E° 
where 

M82 = MSS” + adM..oMSS” — MES”) + bAMSS? — MES” 
with 


(2t—1) 


’ 


M..o = (I — oR,7) [11 — w)I + oRJ(T — oR.) [(1 — w)I + oR,”). 


If the spectral radius of M,... is a[M...] = 1 — 9 (» <1), and if we define 
the rate of convergence of the scheme (9) as 


R(MS2) = —J, in alMS2 
then for large ¢ and optimum a, and b, , we have 


(11) R(MS2) = Vo. 


In what follows we will first establish the convergence of iteration scheme (9) 
and then investigate the dependence of 7 on ¢« for some typical reactor diffusion 
problems. 
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3. Convergence Theorem. A minor manipulation shows that the successive 
overrelaxation operator in equation (5) can be expressed as 


(12) P."LewP = I — wo(I — wR,)*(1 — B). 
Since the spectral radius of B is less than unity, (J — B) is positive definite and 
hence has a positive definite square root. We shall continue to use the notation 


f(A) to denote the spectral radius of A and we shall use || A ||, to denote the spec- 


tral norm of A defined by || A ||, = +/f(AA7) for real A. We now prove 
Lemma. The spectral norm of (J — B)'’P,"L,,.P.(1 — B)~” is less than unity 
for 0 < w < 2. 


Proof. By definition the spectral norm of (J — B)'”’P,"L,.. P.«(I — B)’ 
is the square root of the largest eigenvalue of 


Noo = (I — oT — B)'?(1 — oR.) (U1 — B)"”] 
‘(I — oJ — By? -~ eh") - B)"), 


Obviously, the eigenvalues of N,,,. are non-negative. What has to be shown is that 
they are less than unity for 0 < w < 2. But 


New = I — o(T — B)'?( — oR.) 
(I — oh,”) + (I — oR,) — wo(1 — B)\U — oR,7) “(1 — B)"” 
= I — o(2 — w) (J — B)'*( — oR.) “(I — oR,") “(I — BY” 
= I — w(2 — w)((J — B)'*( — oR) “(I — B)'* (1 — oR)’. 
For 0 < w < 2 we see that we can write 
New = I — PewPs 
where P,,. = ~/w(2 — w)(I — B)'?(I — oR,)™. 


Since the eigenvalues of a real nonsingular matrix times its transpose are real 
and positive, we see that the eigenvalues of N,,.. are less than unity. Hence the 
lemma must follow. 


TuHeEorEM*. The spectral radius of the product of a finite number of successive over- 
relaxation operators || P."Le,uP., is less than unity for 0 < w < 2. 


Proof. The proof is obvious since 


all] P.7L...P.} = a{(1 — B)'[[] P.7L...P.|(1 — B)*"} 


a{] [Ud — B)'’P."L,,.P.(1 — B)~*} 


and hence 
all] P."L...P.} < I] || 7 — B)'?P."L.,.P.(I — B)*”? ||, < 1. 


by the Lemma. 


CoroLuary. The eigenvalues of the symmetric successive overrelaxation operator, 
M,.«, are real, nonnegative, and less than unity for 0 < w < 2. 





* It was pointed out by the referee that this theorem may be derived as a special case of a 
result. by Ostrowski [5]. The proof, however, is different. 
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Proof. M,,.. is the product of two successive overrelaxation operators 
Mee Pi" Lj oP sP.' Le oP 


where the ¢ ordering is the converse of the o ordering. Hence from the theorem, 
the spectral radius of M,,.. is less than unity. That the eigenvalues of M,,.. are 
real and positive follows from (using equation (12) ): 


(I — B)'?M,..(1 — B)? = {1 — o( — B)"*U — oR,") "U1 — B)" 
“(I — o(T — BY’? — oR)“ — By"), 


and noting that this operator is the product of an operator times its transpose. 
Thus we have shown that polynomial extrapolation is applicable (since the 
eigenvalues of M,,, are real) and that the iteration scheme (9) will converge. 


4. The Optimum Extrapolation Parameter. We define the “optimum” extrap- 
olation parameter w for the symmetric successive overrelaxation scheme for a 
given ordering o, as that value of w for which the largest eigenvalue of M,.,. is 
minimized. From equation (12) we see that 


M.o = (I — (I — wh,") “(1 — B) I — o(I — oR.) “(1 — B)| 
I — (I — oR,”) {U1 — oR.) + (1 — oR”) 

— w(I — B)\(1 — oR,) “(1 — B) 
= I — w(2 — w)(1 — oh,”) "(I — wR.) “(1 — B). 


(13) 


Thus we see that the optimum w is that value of w for which the largest eigenvalue of 


! — oB + wR, R, 


(14) w(2 — w) 





is minimized. We note that if X,.° is the largest eigenvalue of (14) and @,° the largest 


eigenvalue of M,... 
0 1 
a 


Moreover, the eigenfunction y,” corresponding to i,” is the eigenfunction of M,.. 
corresponding to @,.". Let 


2 T 
aa-looB+ ok Re —— 
(15) w(2 — w) 


ow =wtio; Yo =’ + Ave. 





We then have 
A(w)ha’ = ru Cy.” 
A(w')Yo+ = AU -CYe-. 


If we take the inner product of the first equation with y¥&- and the second with 
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¥. and subtract the two resulting equations, we obtain, after some manipulation, 


(ur — rw) (Wu Che) = + (ve, — v.') bw 
(16) 





2 2 
+ (ve, 54 ve — (Ay. ,[A — r.’ClAy.) + (d0*) 
1G) 


where the definitions of @A/dw and 3°A/dw’ are obvious. From equation (16) we 
see that ,° has a stationary value for that w for which (y,", (0A/dw)y,") = 0 
is satisfied. This is equivalent to 


2 
(17 = ——— ; P, = 1 — 2t, + 4k 
) eke to + 


where 
To = (yo, By.’); k. = (Yo, R.R."¥.). 


It should be noted that P, > 0, since (y, (J — 2R,|[7 — 2R,."W) > 0 unless 
y =0. Thus 0 < w < 2, and convergence is guaranteed. Moreover, when w 
satisfies (17), a calculation shows that 


aA 
(ve, wy v.') > 0. 


Since A and C are symmetric, A — \,"C is negative semidefinite, hence the third 
term on the right-hand side of (16) is nonnegative. Therefore, when relation (17) 
is satisfied, d,,’ is minimized. For the value of w satisfying relation (17), the spectral 
radius of the symmetric successive overrelaxation method is given by 


a 


om _ w(2 — w)(1 — 14) ta VP. 
(18) 0 * l r 1 — WT w + wk, x 1 + 1 — To 
VP. 


A procedure for estimating the optimum w might be to use a trial function in 
relation (17). In obtaining the few numerical results presented later, this procedure 
seemed fairly insensitive to the trial function. 








5. Comparison with Successive Overrelaxation. For w satisfying equation (17), 
and with 1 — r, = 6, < ~/P, we obtain for 7 in equation (11): 


26 
(19) 7 = ——. 
VP. 


Since k, S land1l — 7, = « = 1 — @(B), where g@(B) is the spectral radius of B, 
we see that the right-hand member of equation (19) has a lower bound 


2e 
20 .» =—. 





0 
— 








rr 
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So we see that for some problems symmetric successive overrelaxation may have a 
convergence rate satisfying 


)\~ 2 
(21) R(MS2) = Ae :. 


In comparing this equation with equation (8), we see that it may be possible for 
successive overrelaxation to be about three times as rapid in convergence as the 
symmetric successive overrelaxation. The lower bound in equation (2) was ob- 
tained when we placed k, = 1. If we assume that P, = O(e) and 6, = O(¢«), 
then 7. = O(¢”); that is, the symmetric successive overrelaxation method will 
converge much faster than the successive overrelaxation method. However, for this 
particular case, we see that we need 


(22) ke = } + O(c). 


6. Numerical Results. An experimental program was written by Miss B. D. 
Baldwin for the IBM 704 to compare the symmetric successive overrelaxation 
method with the successive overrelaxation method in numerically solving the two- 
dimensional diffusion equation 


(23) —VDV¢ + Ag = S; A20; D>0 


in a finite region and with appropriate boundary conditions. 

The Jacobi spectral radius @(B) and fundamental eigenfunction were determined 
by iteration. The successive overrelaxation convergence rate was calculated from 
g(B). The fundamental eigenfunctions of the symmetric successive overrelaxation 
method for various values of w were determined by iteration. For each w the re- 
suling eigenfunction was used as a trial function in equation (17) to calculate 
wopr - Problems representative of typical reactor diffusion calculations were ex- 
amined. The following results were typical: 








: me e : Symmetric Successive Overrelaxation with 
Successive Overrelaxation Polynomial Extrapolation 








| | | = | Sweep 
| 


Ratio 
iB) wort | Rien) Lae ke wort | ROMS) | “ 
RMS?) 
Cale.1.......| 0.988 | 1.73 | 0.27 | 0.988 | 0.334! 1.25 | 0.20 1.35 


Cale.2....... 0.995 | 1.82 | 0.18 | 0.995 | 0.313 1.32 | 0.10; 1.80 








The last column of the table gives the estimated ratio of convergence rates of sym- 
metric successive overrelaxation to successive overrelaxation sweeps. 

In general then, Sheldon’s results do not follow. However, it can be seen that in 
special cases his results will be valid. If we take the coefficients D and A to be con- 
stant and if we use equal mesh spacing in numerically approximating the differen- 
tial equation (23) then it can be shown that equation (22) will be valid. 
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An Improved Eigenvalue Corrector Formula 
for Solving the Schrodinger Equation for 
Central Fields 


By J. W. Cooley 


1. Introduction. The wave equation for the nuclear motion of a diatomic mole- 
cule, in the Born-Oppenheimer approximation, is one which is encountered fre- 
quently in quantum-theoretical calculations. Numerical methods for its solution 
have been developed and used [1, 2, 3, 4] over many years for atomic problems 
where the potential is one obtained by Hartree-Fock self-consistent fields or the 
Thomas-Fermi-Dirac statistical field methods. Only relatively recently have com- 
putational techniques and the application of electronic computers enabled one to 
obtain accurate theoretical internuclear potentials at enough internuclear distances 
to calculate the wave functions for the motion of the nuclei and use them to ob- 
tain averages, over the nuclear motion, of molecular properties. 

The present investigation is concerned with obtaining an accurate method for 
calculating the nuclear wave functions and vibrational-rotational energies of 
diatomic molecules with some economy in the number of values of the internuclear 
potential required. An improved formula for the correction of trial eigenvalues, 
which does not depend so much for its accuracy upon 1. smallness of the step- 
size in the radial coordinate, and an analysis of the conve,_ nee of the procedure 
are given. A computer subroutine was written and numerical results obtained from 
it are described for a case where exact analytical solutions are known. 

In what follows, the vibrational quantum number v, v = 0, 1, 2, --- , will be 
used as a subscript to index the eigenvalues Z, with the usual convention that 
Esk sks-::-. 

2. Method of Integration. The Schrédinger wave equation for the motion of the 
nuclei of a diatomic molecule, regarded as a symmetric top, can be expressed in 
polar spherical coordinates R, 8, & of one nucleus relative to the other and one 
additional coordinate ¢ giving the angular orientation of the electronic charge 
cloud about the internuclear axis. It has been shown [5] that the wave equation 
can be separated into its angular and radial parts and its solution may be expressed 
in the form 


Vv = R* P(R) Yau(®, ®, o) 
where P(?) is the solution of the one-dimensional radial Schroedinger equation 
P®(R) = (U(R) — E) P(R), 
P™ = d"P/dR", 


(2.1) 


the Y,,’s are the hypergeometric functions, A is the quantum number for the 
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z-component of electronic angular momentum, and the radial internuclear potential, 
U(R), is of the form 


U(R) = [J(J +1) — AR? + Z, 2 R* + E.(R). 


The second term is the electrostatic Coulomb repulsion energy of the nuclei and 
E.(R) is the electronic energy obtained by solving the electronic wave equation 
for each fixed internuclear distance R. The boundary conditions for (2.1) are 
P(0) = 0, P(R) bounded. 

To get a difference equation whose solutions approximate the solutions of (2.1) 
we let 





R; = th, i= 0,1,2,---,n+1, 
(2.2) P; = P(R), 
U; = U(R,). 
By dropping terms of the fourth order and higher in the series 
20 2n* 
(2.3) Pin + Pin = > cay P™ 


and by using the differential equation to replace P;”, one obtains the simple 
integration formula 


(2.4) Piga + Pia — 2P; = W(U; — E)P; 


which has an error of approximately (h‘/12)P;“. A higher-order integration 
formula, which does not involve any more values of P; , can be obtained by sub- 
tracting h’/12 times the series 





(2) (2) = 2hn™ (2k+2) 
(2.5) Pin + Pia = & (2k)! P; 
from (2.3). Then, dropping sixth and higher-order terms in h gives 
(2.6) Yio + Yiu — 2Y; = W(U; — E)P; 
where 
(2.7) Y; = P; — (W’/12)P.° = [1 — (W/12)(U; — E)|P;. 


The error in (2.6) is approximately —(h°/240)P;. The integration formula 
(2.6), usually attributed to Numerov [6], involves the extra calculation of the Y,’s 
but reduces the-number of points where values of U(R) are required. 

In the range E > U(~), corresponding to unbound states of the two nuclei, 
solutions of (2.1) exist for all E and can be approximated by simply using (2.4) 
or (2.6) to integrate outward, starting with the boundary values 


(2.8) P, = 0, P, = a small arbitrary number. 


For E < U (@), the two nuclei are bound together and solutions of (2.1) exist 
only for a set of discrete values of E, the eigenvalues of the problem. The method 
for calculating these eigenvalues and corresponding solutions, or eigenfunctions, 
is the subject of the present work. 
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The boundary condition, P(R) bounded, is approximated, as usual, by the 
conditions 


P41, = a small arbitrary number 
P, = Pnat exp (Rass Ona — E—R, VU, — E). 
The second of these conditions results from the assumption that, at R, , U(R) is 
slowly approaching a constant. 
The usual numerical procedure is to provide a first estimate of E and integrate 


outward from R = 0 to some point R#,, , using starting values (2.8) and integration 
formula (2.4) or (2.6). Since (2.1) is homogeneous in P(R), the resulting P; values 


(2.9) 


may be replaced by P;°"* = P,/P,,, i = 1, 2, --- , m. With starting values (2.9), 
the same procedure is used to integrate inward from R,4, to R,, and yields values 
my i = n+ 1,0, ---,msuch that P,” = P,,°“ = 1. Then, a correction to E is 


determined by the difference between the slopes of the two curves at the crossing- 
point R,, and the process is repeated until the two curves meet with the same 


derivative. The correction D(£) [1] (See equation 6, p. 86) is usually calculated by 
using the formula 


(2.10) D(E) = (Pou — Pia) / [ (P(R)}* aR 


where the terms in the numerator are the derivatives at R,, of the curves resulting 
from the outward and inward integration respectively. 

Equation (2.10) has been derived from the differential equation and is claimed 
to give first-order convergence in the error. This means that if a trial solution 
satisfies the differential equation (2.1) at all points except R,, , then (2.10) deviates 
from the true correction by an amount which is proportional to the correction. In 
the next section, it will be shown that the eigenvalues of the difference equations, 
(2.4) or (2.6), are the zeros of appropriate functions F(£) of the trial eigenvalue E 
and that these zeros can conveniently be calculated by the Newton-Raphson 
method. The correction is 


(2.11) D(E) = —F(E)/F’'(E) 

and, for E near an eigenvalue Ep , 

(2.12) D(E) = AE + (AE)’Q + high-order terms in AE 
where AE = Ey — E and 

(2.13) Q = F"(Ey)/2F'( Eo). 


Thus, the convergence of the present method is of second order in AE. In Section 4, 
an analysis of the behavior of D(E) over the whole range of E is given. 


3. Correction formula. The integration formulas (2.4) and (2.6) can be written 
as systems of equations in the P,’s, and Y,’s respectively, as follows: 


(2h* + U,; — E)P,; — h'P, = 0 
(3.1) —h*P,., + (28° + U; — E)P; —hPing =0 
—hP,.4 + (hk? + U, — E)P, =0 
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and 


(2h? + (U, — E){l — (W/12)(U, — E)|¥, —h°Y, =0 


(3.2) —h?Yi4 + (2h + (U; — E){l — (W’/12)(U; — BE) VY: — h°Yign = 0 
—h?Y,41+{ h?>+(U, — E)[l — (h’/12)(Un — E)J VY 2 = 0 
where i = 2, 3,--- , — 1 in both cases. A term involving EF has been omitted 


from the nth equation in both (3.1) and (3.2). It can easily be shown that if one 
has, as one should, started the inward integration at a large enough R, to justify 
the assumption used in forming the boundary condition (2.9), this omission is 
justified. 

The procedure used here to derive and analyze the correction formula is an 
adaptation of a very general and effective technique, due to Léwdin [7], for caleu- 
lating solutions of the Schroedinger equation. To apply Léwdin’s method, consider 
the vector-matrix formulation of equations (3.1) or (3.2), 


(3.3) MC = 0, 


where 0 is a null vector, C is a vector containing the P,’s or the Y;’s, as the case 
may be, and M is the symmetric matrix of coefficients which, for the present, may 
simply be regarded as functions of EZ. After an outward and an inward integration 
with a trial value of £Z, all equations except the mth are satisfied. Now, assume 
that the first row of M and the first element of C correspond to the mth equation 
and the mth variable, respectively, in equations (3.1) or (3.2). Then, by parti- 
tioning the first index of M and C, the result of the integration can be expressed, 


be a ( 1 ) tong 
(3.4) ye 
Ma M.. ie 0 


where F(E) is the amount by which the mth equation of (3.1) or (3.2) is not 
satisfied when integrating with the trial value Z. It is assumed that the first ele- 
ment of C, which is P,, or Y,, , is nonzero and, for convenience, it is set equal to 1. 
Equation (3.4) may be written in the form 


(3.5) F(E) = My + Mi, C, 
(3.6) 0 = Ma + M..C,. 


The function F(£), whose zeros are the eigenvalues of the difference equations, 
is defined by (3.5) in terms of C, which, in turn, is defined as a solution of (3.6). 
Furthermore, if E is such that | M.. | # 0, then the solution C, of (3.6) is unique 
and F(£) is uniquely defined. An expression for F’(E) is obtained by differentiat- 
ing (3.5) and (3.6), with respect to EZ, 


(3.7) F'(E) = My + MC. + Mi, C.’ 
Mi, + M.. Cc. + Maa pay 


(3.8) 0 
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Multiplying (3.8) on the left by C,' gives 
(3.9) 0 = C,' Mi, + C.' Mi. C. + C,' M.. C,’. 


From (3.6) and the fact that M., is symmetric it is seen that if (3.7) and (3.9) are 
added, the terms containing C,’ cancel and the result is 


(3.10) F'(E) = Mi + a Cc. + c.! M:, + c.! M.. Cc. 
which can be written 
(3.11) P’(E) = C'M’C. 


In the two cases considered here, only diagonal elements of M depend upon E 
so that the correction formula (2.11) can be written 


(3.12) D(E) = —F(E) > CM',. 


In the case of (3.1), the derivatives of the diagonal elements of M are all equal to 
—1 so (3.12) becomes 


(3.13) DCE) = ((—Pnaa + 2Pm — Pmi)h” + (Un — E)Pul / > P?. 
i=1 

It is of some interest to compare this with the correction formula (2.10). If, in 
(2.10), central differences at R,, are used to estimate the derivatives and the 
trapezoidal rule is used to estimate the integral in the denominator, then (2.10) 
is the same as (3.13). Therefore, use of the integration formula (3.1) and correction 
formula (2.10) as described above yields a second-order process in AE. 

In the case of the Numerov integration formula, where equations (3.2) are 
solved, the derivatives of the diagonal elements of M are 


(3.14) Mi; = —[1 — (h/12)(U; — EB)y”, i= 1,2,---,n. 


By substituting (2.7) the correction formula (3.12) may be written 


(3.15) DCE) = ((—Ynu + 2¥m — Yuu)h™ + (Un — E)Pal 2, Pe. 

i=l 
Unlike the case with the simpler integration formula (3.1), this is not the correction 
one would obtain from (2.10) by using the most natural higher-order estimates 
of the quantities in (2.10). 


4. Convergence. When the Numerov integration formula (3.2) is used, a 
rigorous treatment of the convergence of the method would be quite complicated. 
Furthermore, since the difference in the two methods considered is a term of order 
h* in the solution P(R) and since both are second-order processes in AE, one is 
justified in assuming that the essential features of the convergence of the two pro- 
cedures will be about the same. Therefore, the convergence of the method using 
equations (3.1) will be considered here. 
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To get expressions for F(£), F’(E) and D(E), let 


Ja 








Gn! K IK 
a a  e 
iG. K ! 
! K ° ! 
: 9 4K 1 
(4.1) H=| K | KZ Gi .' 
—_-——-:| —— ee ee + 1 — ee Se 
K | | Gnoi K 
et oe 
i | G ee 
; ren ( 
| | K G, 
where 
K = -—h” 
G; = 2h’? + Ui, i=1,2,---,n-—1 
G, =h’ + Un 
and where the blank spaces denote zero elements. Then, for equations (3.1), 
(4.2) M=H-IE 


where I is the identity matrix. The symmetric matrix M,, can be diagonalized by 


an orthonormal matrix U, whose columns are eigenvectors of both M,. and Ha. . 
Thus, 


(4.3) M.. = UAU' 


where A is a diagonal matrix with the eigenvalues of M,,. on its diagonal. Then, if 
no eigenvalue of M,, is zero, a solution of (3.6) exists and can be written 


(4.4) C, = —M2. Ma = —UA™'U'M,: . 
Letting 
(4.5) V = U0'M., 


and substituting (4.4) for C, in (3.5), one obtains 
F(E) = My — V'a'V 
= My —_ y af ra. 


(4.6) 


As one can see from (4.2), the diagonal elements of A are 
(4.7) wy =e—E 


where e,,v = 1,2,---,nm — 1, & < €y4;, are the eigenvalues of H.. , the matrix 
obtained by deleting the first row and column of H. Hence, 


(4.8) P(E) = 2h° + U,—E— >> Ve(e, — E)™ 
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and 


(4.9) F(E) = -1— > V-(e, — E)”. 


The set of points R,, --- , Ra+, Russi, --- , Ra, referred to by the subscript 
“q” in the vectors and matrices above, may be partitioned into the set of points 
R,,--- , Ru, used in the outward integration, and the set of points R,.4;, --- , 
R, , used in the inward integration. Letting “out” and “in” subscripts denote the 
sub-vectors and sub-matrices resulting from this partitioning, H.. can be written 


(4.10) a. = “ , ) 
0 Hin 


where H,.; and H;, are the block matrices lying on the diagonal in (4.1) and en- 
closed by the dashed lines. By partitioning in this manner, (4.3) can be written 


Mout ~~ UourAcutU but 
Mi, _ UinAinUl, 


where the columns of U,.; and U;, are the eigenvectors of H,,, and H;, respectively. 
Therefore, e, ,v = 1,2, --- ,n — 1, are the eigenvalues of H,., and H;, . From the 
form of H, it is evident that the e,’s are the values of E which, on the outward 
and inward integrations, respectively, lead to P,, = 0, which is contrary to the 
assumption that the firstyelement of C is nonzero. 

To show that F(£) is undefined at and only at the eigenvalues ¢, it is necessary 
to show that no V, is zero in (4.8). This is evident since (4.5) and the form of 
M.; imply that the elements of V are —h™ times the elements in the last row of 
U.ut and the first row of U;, . If any such element of U..u: or Uj, were zero, then, since 
M..: and M;, are tridiagonal matrices, the column containing that element would 
be zero. This, of course, is not so, since these columns are the eigenvectors of M,... 
and Min . 

It may be of some interest to show how the characteristic polynomial, 


(4.12) P(E) = |M|, 
of (3.1) is related to F(E). From (4.3) and (4.5), one gets 


(4.11) 


_ [Mn v'| 
lvl 


(4.13) P(E) 


which can be expanded in the elements of the first row and column of the determi- 
nant to give 


P(E) = MuJ[» — DV? I] » 
v v vw 0 
= (My — > Vr") I». 
From (4.6) and (4.7), it is seen that 
(4.15) P(E) = F(E)-|] (e — EB). 


(4.14) 
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Equations (4.8) and (4.9) enable one to determine the general behavior of 
F(E), and D(£). They show, first of all, that F(#), F’(£) and D(£) are defined and 
continuous for all E except E = e,,v = 1,2,--- ,n — 1 and that 


F(E) ~2h°+U,—E for | E | large 


F(E) <2h°+U,-—E  forE <e, 
F(E) >2h?+U,—E for E> ens 
(4.16) F(E) ~—V.(e—E)" forE~xe 
F'(E) < -1 for all E 
F'(E) = -1 for | E | large 


F'(E) ~ —V.(e — E)~* for E ~ e . 


Therefore, in each of the intervals into which the points e, divide the E-axis, F(E) 
is a continuous decreasing function which goes from positive to negative values. 
Hence, in each such interval, F(£) must have one and only one zero. Let these 
zeros be denoted by E,, v = 0,1, 2,---,n — 1 with Ay < aq ,e@ < EB, < @yys, 
€n-1 < En. The index v is, therefore, the vibrational quantum number. From 
(4.16), it is apparent that the correction formula (3.13) has the following prop- 
erties: 


D(E) ~E—-e for E ~ e, 

D(E) + 2h° + U, —E for | EF | large 
(4.17) D(E) <2h*7+U,-—E forE<e 
D(E) > 2h*+U,—E for E > ey, 
D(E) SE, — E for E ~ E,. 


The last condition is a general property of the Newton Raphson method (see 
(2.12)). A plot of F(£) and D(E), derived from the above analysis, may be 
expected to appear as shown in Figure 1. 

A serious convergence difficulty is immediately obvious from the first property 
of D(E). This indicates not only that E ~ e, leads to a gross underestimate of the 
correction for such values of Z, but that the smallness of D(E) cannot, by itself, 
be used as a convergence criterion. The condition E ~ e, can easily be detected by 
the existence of a large F’(£) and an increase in D(E) from one iteration to the 
next. Another difficulty which may occur is a jumping from one branch of the 
F(E) curve to another on successive iterations. This can cause one to miss some 
desired eigenvalues and waste computing effort in converging to eigenvalues which 
are not wanted or which have already been computed. Convergence problems, 
therefore, are related to the distribution of the vertical asymptotes at e, , which 
are determined by the crossing point R,, , and to the selection of initial trial values 
of E. 

To investigate the role of R,, , consider the convergence factor (2.13) 


(4.18) g = Ft) 


 2F'(E,) 
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Fig. 1.—Behavior of F(Z) and D(E) curves. 


of the Newton-Raphson method. For E, ~ e,- , 
(4.19) Q = (e — E,)* 


showing that Q is large and convergence is poor if an eigenvalue E, is near a vertical 
asymptote. It was pointed out above that e, is a value of E for which P,, = 0 
on the inward or outward integration. Therefore, E, ~ e,- is the situation where, 
for the correct solution, P,, 0. In other words, convergence is poor if R,, is near 
a node of the desired solution. Some practical means for selecting the crossing-point 
R,, , and initial estimates of EF are given in the next section. 


5. Application. A computer subroutine was written which, when given a numeri- 
cal potential, an initial estimate of FE, and an e, integrates by the Numerov method 
(equation (2.6)) and uses (3.15) to correct R. After D(E) starts decreasing from 
one iteration to the next, the convergence criterion D(#) & « is applied. 

In order to keep R,, as far as possible from a node of the solution and to acquire 
maximum significance where the solution ‘is large, R,, is selected during the in- 
tegration as the point where P; is largest. Since, for the applications considered, 
this occurs at the outermost maximum point of P;, the inward integration is 
performed first and R,, is taken as the point at which P; stops increasing with 
decreasing 7. The description of the behavior of F(E), given in the previous sec- 
tion, is somewhat invalidated since, there, it was assumed that R,, was held fixed 
while varying E. Instead, the F(£) given by the present program will behave like 
F(E) of the previous section in each E-interval where the variation of E does not 
change R,, . 

The data used here to demonstrate the method was obtained for the case where 
U(R) of equation (2.1) is a Morse potential* [8] 





* The parameters used are: a = .711248, R, = 1.9975, D = 188.4355. These were obtained 
by fitting the Morse potential to a computed potential energy curve for H2*. Energies are 
given in units of 2u atomic units, where u is the reduced mass of the nuclei. 
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Fic. 2.—Graph of calculated D(EZ). Crosses on curve denote successive iterates obtained 
with poor starting values Ey’ and 2,’. 






































TABLE | 
Successive Iterates Obtained wiih Poor First Estimates of E, and E, 
i E,®) F( Eo") E,®) F(E:“)) 
—_ | 
1 — 168.800 00 | —1470 — 168.500 00 456.7 
2 —169.450 00 | —727.9 — 166.581 98 | 220.5 
3 —170.700 20 | —350.7 —163.292 37 | 58.21 
4 —172.904 49 | —113.0 — 160.385 38 | 1.654 
5 —176.702 14 —26.24 — 160.283 89 | .0003 
6 | —178.624 20 —2.042 —160.283 69 |  .1550 X 10~ 
7 —178.797 03 — .1786 X 107 
8 —178.798 56 — .5185 X 10+ | 
9 —178.798 57 — .1228 X 10-° | 
TABLE 2 
Dependence of Eigenvalues on n, the Number of Integration Points 
" Ls Bi E: | Bs | Es 
| | 
50 —178.81052| —160.35850) — 143 .02059| —126.83300) —111.81094 
100 —178.79924 — 160. 28784) —142.79397| —126.31918) —110.86348 
150 —178.79866| —160.28428| —142.78276| —126.29441; —110.81921 
200 —178.79857| —160.28368| —142.78090) —126.29031; —110.81191 
Exact | —178.79850) —160.28182) —142.77990| —126.28824) —110.80832 
(5.1) U(R) = D{l — exp (—a(R — R,))f — D. 


Values of the analytic solution are given for comparison. The range 0 < R S 10 
was used in all cases. 

A graph of the calculated D(Z) is given in Figure 2 for Z in the region of the 
two lowest eigenvalues. The effect of the shifting R,, is evident. As one might expect, 
the discontinuities are large in the region of the e,’s, where the two solution curves 
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cross with large slopes. On the other hand, there is a large region about each eigen- 
value E, where such discontinuities are negligible and where D(E£) is almost a 
straight line of slope —1 and is, therefore, a good estimate of the necessary cor- 
rection. 

To show how the procedure converges when a poor first estimate of E is used, 
the program was given a value of E near e, . As one can see in Table 1 where the 
successive trial values are given, and in Figure 2, D(£) increases for a number of 
iterations and then goes down rapidly. Table 2 gives some of the eigenvalues ob- 
tained for n = 50, 100, 150, and 200 points. Exact values of the analytic solution 
are given for comparison. Table 3 contains some of the values of the v = 0 eigen- 
functions which were obtained with 50, 100, and 200 points. Values obtained from 
the exact analytic solution are given in the last column. Figure 3 contains a graph 
of the potential (5.1) and the wave functions obtained for the v = 0, 1 and 2 states. 


6. Conclusions. The above results indicate that fairly good accuracy can be 
achieved rapidly by the above procedure, but that good first trial values are quite 
important, not only in reducing the number of iterations, but in assuring that no 



































TABLE 3 
Dependence of Wave Function P(R), for v = 0; on n, the 
Number of Integration Points 
P(R) 
R : 
n = 50 nm = 100 | n = 200 | Exact 
1.2 .022 5903 | 022 3875 | 022 3754 | .022 3746 
1.6 .484 0974 .485 8844 .485 9864 .485 9927 
2.0 1.312 2858 | 1.310 2630 | 1.310 1471 | 1.310 1405 
2.4 .732 9700 | .734 7376 .734 8413 | .734 8476 
2.8 .126 3448 .126 4907 .126 4998 .126 5004 
3.2 .008 9574 | .008 9544 | .008 9541 | .008 9541 
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Fie. 3.—Plot of potential U(R) and v = 0, 1, 2 solutions. 
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desired solutions are missed. Good first estimates of E, can be obtained by fitting 
an analytic potential for which eigenvalues are known and using it to obtain first 
estimates of the lowest eigenvalues. An extrapolation of calculated eigenvalues 
can then give an approximation to each new eigenvalue. 

If it is inconvenient to obtain first estimates of E in this manner, one can have 
the computer subroutine perform just one iteration for each of a series of values of 
E to get a D(E) curve. Then, in each interval where D(E) changes from positive 
to negative, there is an eigenvalue for which a first estimate can be obtained by 
interpolation of D(£). A count of the nodes in the solution, given by the subrou- 
tine, is a check on the vibrational quantum number assigned and assures that no 
solution is missed. 
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Integration of the General Bivariate Gaussian 
Distribution over an Offset Circle 


By A. R. DiDonato and M. P. Jarnagin 


1. Introduction. An efficient method is described in this paper for the numerical 
evaluation by a high-speed digital computer of the integral of any uncorrelated 
elliptical Gaussian distribution over the area of any arbitrarily centered circle in 
the plane. Actually, by introducing at most four linear transformations of coordi- 
nates, the present method can be used to cover the general problem of integrating 
any correlated or uncorrelated bivariate normal distribution over an arbitrary 
ellipse in the plane. 

The basic integral to be evaluated is expressible in the form 


> 1 x y 
(1) I -scl. i p(r) exp | - (5 + v) dxdy, 


where the exponential represents an uncorrelated bivariate normal distribution 
centered at the origin with standard deviations ¢, and o,. A point target T is 
assumed at some arbitrary point (h, k) in the Oxy plane. The distance r is ordinary 
distance [(z — h)? + (y — k)*]"”, where (zx, y) is an arbitrary burst point in the 
plane. The function p(r) by definition has the value one if r does not exceed a 
given constant R, and p(r) is zero otherwise. This has the effect of limiting the 


field of integration to a circle with center at (h, k) and radius R. Hence the integral 
can be written as 


k++/R2—(z—h)?2 1 x y 
» regia Sere +Hle 
(2) eo h—-r Jk—/R—G—ht PL 2 \e2 " o,? gee: 


The problem was previously studied in [2], [3], and was mentioned in [4]. The 
more general problem of integrating an arbitrary integrable function of two vari- 
ables over the circle has been analyzed in [5], [6] and other papers. 





2. Transformations of the Basic Integral. The efficient numerical integration 
of the iterated integral of equation (2) is not a straightforward procedure for the 
following reasons: 

(a) The behavior of the integrand in the second integration is such that for 
particular values of the input parameters it attains the shape of a spike with vir- 
tually all of its area contained within perhaps 1/100 or 1/1000 of the entire range 
of integration [h — R,h + Rj; 

(b) The slope with respect to x of the integrand in the second integration ap- 
proaches infinity as x approaches h + R. 

The minimum requirements for the ranges of values of the input variables were 

1 


6 


=<6, oxs*<, o<*<20. 
Cy Oz Cy 


IIA 
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These were considered realistic ranges in the light of likely applications, for example, 
to problems of weapons assessment and other problems of operations research. 
The program which was developed, however, is valid for much wider ranges of 
values, as discussed in Section 8. 

If a normalizing substitution 


(3) u=(x—h)/R 


is made in the integral of equation (2), the range of integration on u is from —1 
to 1. If the resulting integral in u is split at zero and the sign of u is changed in 
the integral from —1 to 0, the integration range can be reduced to the interval 
from zero to one. The difficulty of an unbounded slope of the integrand in the 
neighborhood of the upper limit of integration, here u = 1, is resolved by making 
another substitution, 


(4) t=~SJfl—u 


yielding P in the following form: 


R/o: f° ifh-—RA-? ifh+ Ri —?’)f 
lal ved {en(-5[*=- BG = ‘T) * ap(-;[St ’T)} 
{Erf (du) — Erf (dy) }t de 
where 
aE. k+ Rinf2 -—# 
V 20, . 
(6) Erf (x) = Kl [ e” dv. 
a ba RtV2= 8 _ 
01 /20, 


Equation (5) can be considered as a single integral if the Erf (dy) and Erf (do) 
are numerically evaluated by analytic procedures rather than by quadratures. This 
is actually what is done. A table of the functions Erf (x) and Erf’ (x) is assumed 
stored in the computer at suitably small intervals of the argument. The value of 
either of these functions, for any required value of the argument, is rapidly de- 
termined by a table look-up and the use of no more than four additional terms of 
the Taylor expansion about the nearest value for which the functions are stored. 
The truncation error, or remainder term, is then subject to a rigorous predetermined 
bound, as discussed of Section VI of [1]. The efficient methods of Hammer [5], 
Peirce [6], and others, for integration by quadrature over multi-dimensional do- 
mains, were considered at the request of the referee. It appears evident, however, 
that it is more efficient to take advantage of the opportunity to perform one of the 
two integrations analytically, since it can be done much more rapidly and systema- 
ically in this way than by any numerical quadrature method, and the other by 
Gaussian quadrature for single integrals. A small sample of test cases was run and 
this conclusion was in general verified. See Table 1. 








3. Determination of Reduced Intervals of Integration. The difficulty described 
in the preceding section under (a) can be virtually eliminated by changing the 
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limits of integration on ¢ to é and e, , where 0 S e S & S 1, and to a lesser degree 
by changing the limits of the error functions that occur in equation (5) from dy 
and dy to dy and d, , where dy < db S dh S dy. 

For a given «¢ in the range 0 < « < 1 there is a unique rectangle, S, centered 
at the origin, with sides of lengths 2a(«)o, and 2a(«)c, , which are parallel to the 


x- and y-axes, respectively (see Figure 1), such that 


psa Coola B+ Baw 
(7) i 24020, ac, *—aty exp | . os + o,? cated * ” 


Transforming this integral by elementary methods, the following relation between 
« and a can be derived 


(8) Erf (55) =vVi-<. 

All but ¢ of the entire bivariate distribution falls within the rectangle, by equation 
(7). Therefore, in carrying out the integration over the circle, as indicated by equa- 
tion (5), only that part of the interior of the circle which also lies in the interior 
of the rectangle need be considered. If one integrates only over the common part 
of the interior of the circle with that of the rectangle (shaded area in Figure 1), 


denoting the result by P, , an error less than ¢ in magnitude will be committed. 
Thus the values of ¢ can be restricted by the following inequalities 


R—h-— ae, 


R—-—h+ae 
e ee - FR 
R t 


(9) s R 


lA 


This can be seen by noting that the left and right boundaries of the rectangle are 
the lines x = —ao, and x = ao, , and getting the corresponding value of t, using 
equations (3) and (4). But the entire integration interval on ¢ in equation (5) is 
from zero to one. Hence the actual range of integration [ep , e:| will be the common 
part of the interval, specified by (9), with the interval [0, 1]. The effect of this is 


¥ ) 
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that the limits e, and e; are given by 














/R—h+ az. h — ac, 
(10a) en Trdee R <1 
(10b) a =4 1 if = P <0 
4 h — ac, 
> 
(10e) 0 if p= 1 
( /R—h—ae.. h + as, 
eo = { h 
h + doz . 
(10e) | 0 f Rr 2 1 


In order to keep the interval over which the integration is carried out as small 
as possible, if the analogous difference e; — e’ for o, and k is smaller than the 
difference e,; — e based on o, and h, x, oz , and h are interchanged with y, o, , and 
k respectively. This interchange is always possible because of the similar roles 
played in equation (1) by z and y and their associated parameters. 

The integration interval [du , dj of equation (5) can also be reduced in many 
cases by analysis similar to that employed above for [eo , e:]. The new limits will 
be defined as follows: 








(11a) a a dy if du < a 
(11b) ale") ie a, = 280 
/2 vbyadn /2 
. (e) 
(12a) tte > -F 
2 
dy = } ( *) ms 
(12b) — <<" if du < _ ale) 


where 


oc 
2€ T6E- 


Further, since, by equation (8), Erf [a(e*)/+/2] = ~/1 — e* = 1 — 4e, Erf (dy) 
in equation (5) can be replaced by (1 — ¢/4) if (11b) holds, and similarly Erf (do) 
can be replaced by («/4 — 1) if (12b) holds. Moreover, because of the monotonic 
character of dy and do as functions of ¢ (see equations (6)), increasing and 
decreasing respectively on the interval [0, 1], if (11b) holds for a particular value 
of t, it will hold for all subsequent values of ¢ in the integration, and similarly 
for do, and (12b). A detailed analysis is given in Appendix A of [1]. 

Thus the probability P as given by equation (5) is approximated within an 
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arbitrarily chosen « by P; , where 


i R/o, “ _fP 
(13a) P= 24/8 Sex (fi + fo) [Erf (d,) — Erf (dy))tde 


te ill | SR OT) 
fait he Feo(-3[ 
2 1 a+ Bu 7) 
ss ae (  pecomer | 5 
4. Determination of Conditions for P < 5«, P > 1 — 5e. In order to improve 
the efficiency of the machine computing program, conditions were determined 
such that, if the kill probability were less than 5e, a result of zero would be stored 


without further computation, or a result of one would be stored if P were greater 
than 1 — 5e. The details are given in [1]. 





(13b) 


5. Determination of the Gaussian Order of Integration. The method chosen for 
evaluating equation (13a) numerically is that of Gaussian quadrature. The symbol 
O(G) will denote the order of Gaussian multipliers. The orders stored in the NORC 
program are 6, 8, 12, 16, 20, and 24 [7]. For the ranges of values of the input pa- 
rameters for which the program is designed, it is never necessary to use more than 
16(24) integration points over the entire interval of integration in order to obtain 
results correct to three (six) decimal digits. 

The rule for determining O(G) in the program is empirical, but it was developed 
by an extensive program of test cases, with adequate checking of results, as dis- 
cussed in Section 9. 

Two levels of accuracy were considered, three and six decimal digits for the 
computed probabilities. For each accuracy level, O(G) depends on the magnitude 
of a positive quantity N given by the empirical equation 





3 i R/ey 
(14) N= (e €) [017 Ries + gels | 


In the program for six-digit accuracy (« = 107’), O(G) is given by the following 
rules: , 


— 


f N = 2.75, 0(G) 
if 0.8 < N < 14, 0(G) 


24; if 14 < N < 2.75, O(G) = 20; 
16; if 0.4 < N < 0.8, 0(G) = 12; 
if 0.15 < N < 04, 0(G) = 8; if N < 0.15, O(G) = 6. 


In the program for three-digit accuracy (« = 10“): 


If N = 2.0, O(G) = 16; if 0.8 = N < 2.0,0(G) = 12; 
if0.5 < N < 0.8, 0(G) = 8;if N < 0.5, O(G) = 6. 


6. Computation of the Error Function (Erf (x)) and Its Derivative. The required 
values of Erf (x) and its derivative, (2/+/x)e~, are obtained with a subroutine 
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based on Taylor series expansions [8] about stored values of these functions, and 
taking account of theoretical bounds on the error terms in the series. The numerical 
values of the error function and its derivative are stored at intervals of 0.01 from 
x = 0 to x = 3.85. Therefore this subroutine is very fast in that it requires only 
two terms of the Taylor series for four-digit accuracy in the computed functions, 
and only four terms for nine-digit accuracy, including in each case the rapidly 
accessible stored value. These accuracy levels are acceptable in the programs for 
three and six decimal digit computed probabilities, respectively, as discussed in 
detail in Section VII of [1]. 


7. Error Analysis. A detailed error analysis of the computations is omitted 
here. However, it is shown in Section VII of [1] that for every case within the per- 
missible ranges of the input parameters, the total error in the computed prob- 
ability, due to round-off, truncation of series, and all other causes, is less than 5e. 
The quantity «¢ is taken as 10“ and 10~’ in the programs for three and six decimal 
digit accuracy respectively. 


8. Nature and Scope of the Computing Program. The program is constructed 
so as to accept the five input parameters R, oz, o, , h, k, even though there exist 
only four independent combinations. The limitations in the permissible ranges 
for the present program are 


IA 


2/Gy S 15; 
600; 
600; 


IIA 
lA 


h/o-z 
k/ Cy 
R> 


1s 
0 

(15) 
0 


IIA 
IIA 


S 


These ranges greatly exceed the originally contemplated requirements (Section 2 
above). The primary reasons for this extension were to analyze the function P 
thoroughly and to allow for possible extensions of the realistic ranges in the future. 


9. Discussion of Results. Since the order O(G) of Gaussian multipliers for a 
given case is determined by an empirical rule (Section 5), it was necessary to check 
the validity of the rule with an extensive test program of cases spanning the space 
of permissible values of the input parameters. More than 7,000 cases were com- 
puted in this test program, and the accuracy of the results was verified by doubling 
the number of integration points, and also in borderline cases by checking against 
values computed by an independent iterative Simpson’s rule routine. In addition, 
the tables of values given in [2] and [3] were computed independently by the present 
program. The results agreed, within one unit in the last place retained by the 
respective authors, except for one discrepancy in the case of each author, where 
the difference was more than one unit. In these two cases, as well as in all other 
cases of the total of more than 7,000, it was established that the results given by 
the present program were correct to the specified extent, three or six decimal digits. 

The input values used in the test program are given in Table III of [1]. Briefly, 
o, = 1, o, = 1, 3, 6, 10, 15; h and k range from 0 to 600. For a given set of values 
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of o, , h and k, ten values of R were determined and used in such a way that the 
corresponding values of P ranged from roughly 10“ to 1 — 10°°. 

A table of 10,080 entries is given in Appendix C of [1]. This is an inverse table 
giving R as a function of P, o, (=1), o,, h, k. Six values of P are used, .05, .20, 
50, .70, .90, .95. The values of R were computed such that the specified probabili- 
ties, P, were obtained to within three or four units in the fifth significant digit. 


U. 8. Naval Weapons Laboratory 
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An A Priori Determination of Serial 
Correlation in Computer Generated 
Random Numbers 


By Martin Greenberger 


1. Background. Ever since John Von Neumann introduced the “‘mid-square” 
method some ten years ago [11], users of “pseudo-random” numbers on modern 
digital computers have been generating their numbers as they were required in 
application, rather than drawing them from a table formed beforehand or obtaining 
them from a special-purpose device built expressly for this service. The advantages 
of Von Neumann’s method were its speed, minimal storage requirements, and 
ability to be restarted from any desired point. 

For some years now the mid-square method has been replaced by other recursive 
methods which have increased further the advantages of this general technique of 
generation. The most popular of these has been the “multiplicative congruential”’ 
method proposed by Lehmer in 1951 [2]. An interesting variation of this method, 
which we may call the “mixed congruential” method, recently has received the 
attention of several authors [1, 9, 10]. We adopt the notation of Coveyou, one of 
these authors. The mixed congruential method then may be written 


(1) Ings = An + (mod P). 


Here the modulus P generally equals one more than the largest (fixed-point) 
integer which the computer can store. The multiplier \ and the addend (or in- 
crement) yu are, to a degree, optional parameters of the generator; both parameters 
are positive integers less than P and relatively prime to P. The starting value x» 
is the first term of the integer sequence {z,: 0 < x, < P(n = 0, 1, 2,---)} 
generated recursively by equation (1). And the numbers {z,/P} are the desired uni- 
formly distributed drawings from the unit interval. For the special case of u = 0, 
equation (1) becomes the familiar multiplicative congruential method. 


2. Selection of Parameters. It is not difficult to find values of the parameters 
\ and yw which produce a maximum period in the {z,} sequence for a given P [3-8, 
10]. When P is a power of 2, for example, any odd » combined with any A = 1 
modulo 4 affords the maximum period, equal to P. For » = 0, the maximum period 
is reduced to P/4 and is attainable again with half the odd A, now A = 3 or 5 modulo 
8 [10]. 

Length of period is one standard that has been used in the selection of param- 
eters; speed of generation is a second [1, 3, 4, 10]. But a great deal of freedom still 
remains, and the question of how best to use this freedom never has been fully 
answered. The most frequent solution has been to make a choice of parameters 
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which appears favorable on intuitive grounds*, and then carefully examine gen- 
erated subsequences by the standard statistical tests for serial correlations, 
runs, interval frequencies, and so on. 

If we had a complete understanding of the relationship between the number theo- 
retic properties of A, u, and P, on the one hand, and the statistical properties of 
the sequence they generate, on the other, the selection problem essentially would 
be solved. The fact is that we are still a considerable distance from having this 
complete understanding. Some recent work [9], however, has provided a start in 
the right general direction. We set out now, first to comment on this work, and 
then to develop some additional results which will help in the choice of generator 
parameters. 


3. Serial Correlation. In his article on ‘Serial Correlation in the Generation 
of Pseudo-random Numbers” [9], Coveyou gives the following approximate formula 
for the serial correlation, p(, , 2n4:), between a number and its immediate suc- 
cessor in a full {z,} sequence generated by equation (1). (The symbol = denotes 
“approximately equal to’’). 


_1 6p Mm 
(2) pP(%n, Tny1) = x )P ( fo *) 
Equation (2) is a good approximation to p(n , Zn4:) for \ small compared to P””. 
For \ on the order of P“? or larger, however, the approximation can fail badly, as 
an example will illustrate. 

For P = 2°, = 2" + 1, and » = 1, the correlation given by equation (2) is a 
negligible 2“. The true correlation, however, is found by direct calculation to be a 
very significant .25, much too large to be acceptable. This is explained by the fact 
that values of \ which are very large relative to P have an effect similar to that of 
values of \ which are very small relative to P. This fact is completely distorted by 
the erroneous implication of equation (2) that the larger \, the smaller the magni- 
tude of the correlation. 

As a consequence, the usefulness of equation (2) as a standard for selecting 
parameter values is limited. In what follows, an exact derivation of p(t, , tn4:) 
will provide us with a correction term for equation (2) which overcomes this limi- 
tation. The derivation begins along the same general lines as Coveyou’s skillful 
approach. 


4. Derivation. We shall assume throughout the derivation that ’ and yu are 
restricted to values for which the period of the {z,} sequence is P. The sequence 
therefore will consist of all integers from 0 to P — 1 (chaotically disordered), with 
each integer appearing exactly once. Let the symbol E designate an expected or 
mean value over the full sequence. Then p(x, , %n4:) may be expressed as follows: 





PBN E(2, 2 Tn 1) sie [E(2,)) 
mm p(n» Batt) — Fat) — ECE 
where 
(4) E(2.) = oa i Aol 





* One line of intuitive reasoning has led to a proposal that \ be on the order of P'/? (3, 4). 
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(5) E(2,2) = 1 > r¥- (P — 1)(2P — 1) 


Pi 6 





To evaluate E(x, , n4:), we equate z,4, with the remainder r, in the following 
equation 


(6) Wn tw = GP +r (n = 0,1, ---) 


where g, and r, < P are non-negative integers, uniquely determined for each z, 
by the Euclidean algorithm. Dispensing with the subscript n in equation (6) for 
convenience, we now may write 


P—1 


v 


1 
E(an » Zn41) = P=, P 


x(Ar + uw — gP) 


In Ta+i = 


IM 


(7) 


~ 


—1 


= E(x”) + wE(zx) — 


IM 


rq. 


Assume for the moment that u 2 X and let z, , or z, in equation (6) take on 
consecutive integral values starting with 0. Then g increases from 0 to A, in in- 
crements of 1, and each gq has associated with it approximately P/A consecutive 
integers x, as well as the same number of integers r. Let #, be the smallest r as- 
sociated with a given g. Then 7 = yu, and for g = 1, #, < is given by 


(8) f, =p— oP (mod i). 
From this it follows that the 7, are distinct and (7; , 72, --- , ®&) is a permutation 
of the integers (0, 1, --- ,A — 1). 
It may be shown that 
SP oP wt nF Mo, +20 
z=0 q 2 2 2r 2 22 - q 2)? q=1 q 
~ 1l< A+ Ae P< 
= 2 
~ ae ee" + >; La- yd Te 
and hence that 
P-—1 2 2 : 
~~ AP AP _u i: wi pP r 
im yu-s- poate Rta tD 
10 
aad Be Pi oe 
2 12d 4 12d » Xe 
where 
X 
(11) S= Dory. 
q=1 


5. An Exact Formula for p. By combining equations (3, 4, 5, 7, and 10), we may 
now write a formula for p(z, , 2n4:). In so doing we make the assumption that P 
is so large that any term whose order of magnitude is 1/P or less is negligible. 
This assumption leads to 
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; 2 
(12) p( In 5 Zn41) 21%, -%)+2(5 -3). 


A AP P 


Since P is generally at least one billion, the approximation (=) is an equality for 
all practical purposes. 

The first two terms of equation (12) are seen to agree exactly with the approxi- 
mation for p, equation (2), given earlier. The final term of equation (12) thus 
constitutes a correction to the earlier result, and it will be interesting to examine 
this term carefully. It is not difficult to show that 


= 3 


r 


Xr 2 2 
. —-1)sSs- -—1 
(13) 5 )eSszQ ) 
and hence, because of the magnitude of P, 
r 12/S r r 
ao «= & aot ae ooo & Be . 
(14) eS ( : *) ae 


Thus, the correction term is confined to an interval extending a distance \/P 
on either side of zero. At the mid-point of the interval, when the correction term 
equals zero, equation (12) reduces to equation (2). Equation (2) also suffices 
when 4 is small relative to P’’, or more precisely, when \/P < 1/X. But for \ on 
the order of P“? or larger, the correction term may predominate and the complete 
equation (12) must be used. 

If we repeat the line of reasoning leading to equation (12) for the case un < X, 
we find that 


; ie 12/S r 
(15) P(2n 5 Tay) = x ss P (5 *) 
where 

ar 
(16) S= Lira 
q=1 


and inequality (14) again obtains. In a strict sense the notation should be modified 
to distinguish the S of equation (16) from the S of equation (11). However, the 
context (i.e., whether » = \ or » < Xd) always identifies which of the definitions 
for S applies, and, in any case, the two definitions cannot give results for p different 
by more than a negligible amount. The implication of equation (15) is that equation 
(12) still applies for the case of » < X. For this case, however, the middle term is 
insignificant and the parameter u appears only insofar as it influences the value of S. 


6. Examples and Applications. To evaluate S for some specific cases, let U be 
the unique non-negative integer less than \ for which 


(17) U=u (mod \) 
and let 
(18) V= : <1 
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In what follows, » is assumed to be a positive odd integer less than P = . 
First, suppose that = 2“ + 1. Then it may be shown* that 


(19) S=(V?—V+4)2™. 


The referee has pointed out that equation (1) is the reciprocity formula for Dede- 
kind sums (after an elementary transformation). He has kindly supplied the 
reference to [12]. For V close to either 0 or 1 (e.g., 1 = 1), equation (19) leads to 
S + 42 and the term involving S in equation (12) dominates. The other terms 
are negligible in comparison and equation (12) reduces to 


(20) p(in, Inxs) = (GE — 1) =F 
which agrees with the result given earlier. 


Next, suppose that \ = 2" + 1. Then, if U < 2” 


, 


(21) S = (—2V? + V + ¥y)2” 
whereas if U > 2”, 
(22) S = (—2V* + 3V — yy)2”. 


Equations (21) and (22) may be used with equation (12) to appraise different 
choices of uw. Let is examine u = 1 as an illustration. Thén U = 1, V = 0,S = 
(1%)2™, and 


(23) p(tn, ny) = 2” + 12(2-")( — 2) = 2" -2” 


from which we infer that p< 2”. 

It is not correct to conclude from this evidence alone, however, that a combi- 
nation in equation (1) of P = 2° = 2" + 1, and w = 1 furnishes an acceptable 
pseudo-random number generator. As a matter of fact, it does not [10]. For one 
thing, the first several hundred numbers generated by this combination, starting 
with xz» = 0, are all less than P/2. 

Now consider \ = 2"’ + 1. With this A, 


(24) S = (-—V* + V + $)2° 

so that for V close to either 0 or 1, S = (%)2” and 

(25) (Xn, Xngi) = 3(w +2” — p-2” + 1)2™. 

Thus, when » = 1, p = 3(2””); and when uy simultaneously satisfies » = (1 + 


2°“?)9" and U = 0,p x2”. 





* In each of the examples discussed, special methods have been used to evaluate S, and 
the details are available upon request. The methods differ from each other somewhat, de- 
pending upon the values assumed by the parameters. It is actually possible to develop one 
general technique, a reciprocity type of reduction method, for solving all of the examples. 
S may be expressed in terms of an S’, with \ = P’, by retracing many of the same steps which 
led to equation (12); and so on. The general technique has not yet proven as advantageous as 
the special methods in particular applications, but is of theoretical interest nonetheless. An 
elegant formulation of the general technique recently was developed by Coveyou in a private 
memorandum to the author. 
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Finally, consider ) sufficiently small relative to P“? so that \/P « 1/d. Then 
inequality (14) indicates that the correction term in equation (12) is insignificant 
compared to 1/A and therefore may be dropped. ig ane 


1 (6y" 
(26) plas aan) + + (98 — S# 4-1). 
It follows that p = 1/ when u « P, whereas p < 1/A when » = P(1 + 3°%”)/2. 
Here \ and yu are restricted to values which afford the full period. One such selec- 
tion of P = 2”, namely \ = 2’ + 1 and » = 1, has been tested empirically and 
proposed as a suitable generator [1]. 


7. Higher-Order Correlations. As is pointed out by Coveyou [9], equation (12) 
ean be adapted to give correlations of the {z,} sequence with lags greater than 1. 
To accomplish this, we define non-negative integers A,, and yu, , both less than P, 
by the congruences 


(27) Am =X" and pm = ora (mod P) 
where m = 1,2, --- . By applying equation (1) repeatedly, we may write 
(28) Lnam = Amin + Um (mod P) 


which is identical to equation (1) in the special case m = 1. 

Equation (28) generates a sequence in which 2,4 is the immediate successor 
to z, . If \», and u,, are such as to make the period of this sequence P, then equation 
(12) with » replaced by X,, and uw replaced by yu, may be used to evaluate p(z, , 
In4m)- This is then the correlation with lag m of the original {z,} sequence. 

Thus, the evaluation of mth-order correlations introduces no new problems 
so long as the selection of \ and yu leads to pairs (Am , um) Which produce full periods. 
In this regard we note that with P a power of 2, » odd, and A = 1 modulo 4,things 
work out well. For then u,, is odd and X,, = 1 modulo 4 for every positive integer 
m. Thus, the sequence generated by equation (28) does have period P [10], and 
equation (12) may be used to calculate a serial correlation of any order. 
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Expansions of Hypergeometric Functions in 
Hypergeometric Functions 


By Jerry L. Fields and Jet Wimp 


Abstract. In [1] Luke gave an expansion of the confluent hypergeometric func- 
tion in terms of the modified Bessel functions J,(z). The existence of other, similar 
expansions implied that more general expansions might exist. Such was the case. 
Here multiplication type expansions of low-order hypergeometric functions in 


terms of other hypergeometric functions are generalized by Laplace transform 
techniques. 


1. General Expansions. The generalized hypergeometric function ,F,(z), [2], 
is defined by 


} = y %, te » Ap ———— . - 
_— ge ? kl’ 
(1.1) (b;), 


Tr 
where (oc), = Me +») 
| T'(o) 
We assume that no a; is equal to any b; and that no b; is a negative integer. For 
ease in writing, we employ the contracted notation 


4 ay| \_ (a), 2 

(1.2) Fads) = Pa (%|s) =F Go 5 

Thus (a,); is to be interpreted as bE (a;), and similarly for (b,), . Considered 
as & power series in z, »/’,(z) has a radius of convergence equal to infinity if p < q 
and equal to unity if p = q + 1. In general, ,F,(z) is not defined if p = q + 2. 
However, in this paper, we shall say that ,/,(zw) is equal to another series if the 
coefficients of (w)* on both sides are equal, regardless of the relationship between 
p and q. 

Our first expansion is 


ap, c, | ine —_ (ay) n(a@)n(B)n(—z)” 
reP ee (5rrg 20) = & (ba)aly + n)an! 


i i (" +a,n+6,n+a, :) di ic n+ ¥,¢ v). 





(1.3) 


2n+7y+1,n+), a, B, d, 
From (1.3) we may derive other expansions by confluence: 


Gy, C,| _ © (ap) n(ae)n(—z)” 
resPase (5076 |e) = 3; Hedaledel—2)" 


n+an+a —n,n + Y, Cr 
F ; 4 z) oF, , mS ee 
R eee bo +y+ijn+b,|\7)™ a, d, , 
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oul ete (57? ¢| 20) = , (ay) n(a)n( —z) oul, (" +a,n +a, :) 
q9“s 





n—0 (bg)n 1! n+ b, 


XK nF ogs eee | w) , 
Gp, Cr| a = (ay)n(—z)" (" + a,| ) 
(1.6) ones (; ,d, ew) ef p> (b,),n! Fe n+ b, f 


x ral, er ” w) - 


For example, if in (1.3) we replace r by r + 1, z by 2/c,4;, set 8 = c,4; and 
let 8 — «, (1.4) results; (1.5) and (1.6) may be similarly derived; (1.6) is a result 
given previously by Meijer, [7], 1953, page 355, but it is the only result above which 
can be deduced from his work. 

Equation (1.3) is proved by induction on p, q, r, and s. The case p = gq = r = 
s = 0 reduces to 


(1.5) 


x 


my) zw _ > (a)n(B)n(—z)” 1 (" + a,n + B! Y n,n+¥ ) 
(1.7) e€ -2-cta * n+y+1 zZ oF a, B w , 


a result given by Luke, see (1.8) of [3]. The proof of (1.3) then proceeds by Laplace 
and inverse-Laplace transform techniques. Assume (1.3) true for p, q, r and s, 
multiply both sides of (1.3) by (w)*", take the Laplace transform of both sides 
with respect to w, and we obtain, see (17), page 219 of [4], 


7 —hw_ o—l 7 ap, C, = T'(c) Gp, Cryo \z 
[ ew ptrF ots o id, zw) dw = Ay" ptrtil ote ( b,, ds ‘) 


me I'(c) = (ay) n(a)n(B)n( —z)” 
18) = Ay 237 Goaly + n)onl 





F n+a,n + B,n + Gp), 
pre et"! On + 74+ 1,2 4+ b, 


_ —n,n + ¥,¢-,¢)1 
x naFo( a, B, d, 5): 


The induction on r is completed by replacing 1/A by w in (1.10). The induction 
with respect to s is effected by multiplying both sides of (1.3) by (w)’, letting 
w = 1/\ and applying the inverse-Laplace transform, see (1), page 297 of [4]. 
If the above Laplace transform techniques are applied to z instead of w, the in- 
ductions on p and q can be similarly effected. 


2. Specialized Expansions. In this section we give several interesting cases of 
(1.3)-(1.6). 
From (1.3) we have 


¢, | . 2” = (y)n(—z)” 
si (j ew) - asa -oy Sans a— a 


—n, Nn + 7; Cr , 
x rsaP os (2 (y + 1)/2, d, i w) , 








(2.1) 


where we have used the relation, see [2], page 101, (6), 
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(2.2) [s+a- aa) =F, (°52*|:). 


The expansion 


Cr 
*F, (5 











e!? 9? 
si) = {BoA Io(e/2) 


+2 > (-1)” i (n) ey rol o4t ce aot oy, “| w) Tete (5)} 


(2.3) 


n=l 4 + 7’ d, 2 
follows from (1.4), where J,(z) is the modified Bessel function of the first kind. 
Here, use is made of the expression 


(2.4) Ty +1) (2 = (2/2)’e™* iF; Elaes | 22). 


Also from (1.4) we get 


a rt 6+ 1,4, | 
Fa (S| 20) = aon (, +6 +2,b,| :) 
(2.5) Ed (ap) n(z)” Pp ( n+6+1,n+a, |.) 
mi (bn tatBpBt+i,” “\Qn+a+684+2,n+5, 
x P, "(20 — 1), 








where 





(26) P,"(2) = (-1)" At Be af. (— Be hey At res) 


The P,‘*” (x) are known as the Jacobi polynomials, see [8], and reduce to the 
Chebyshev polynomials of the first kind if a = 8 = —}. 
Equation (1.5) yields the following expansions: 


By -- —N, C,| 
(2.7) F431 Ay ;|=) =e D Se - nan ( na 'w); 
| 


(28) .F, (5 |») = (i= ay Se(-4) waPon (Zo w); and 


¢ |  \N  & Sr(b+n,2) —n, C,| 
(2.9) Fa * + 3| nu) _ (z)8 ee poe nF ( 5, d, ra ? 


where 





ll 


l(a, z) = T(a) — ft et dt, 


(2.11) RAtd zw) = > (ae)al me)" or, (" aieeindny ‘ajind *|2) L,*(w), 


(2.10) 


ll 
a 
ie 
> 
Y every, 
2! 
+|= 
—_— 
| 
nN 
‘ee 





n=0 eee n as b, 
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where 
a oat (a + 1), —n bas 
(2.12) L, (zx) _ ~ ane T iF; 1+ a|” 


The L,*(x) are known as the generalized Laguerre polynomials, see {8}. 

The result (2.7) can also be deduced from the work of Rainville [5], page 267, 
(25), and (2.8) is a generalization of a result given by Chaundy [6] and Meijc. 
[7], 1952, page 483. 


If we use the formula 


(2.13) P, ( + i=) = (> +1) (2) L(2), 


(1.6) yields the expansion 


C, | ow 
Fas (.. 1+] *) 


(2.14) ($y 
= T(1 +8) 2\' > 2 Tasa(2) ea? —nN, C, | 
_ } 2 ——— b+n\Z) roils d \w)}. 


n=0 2 | 


Since the merits of modified Bessel function expansions are well known, we 
consider the convergence properties of (2.3) in more detail. They are determined 
principally by the asymptotic properties of the polynomials 


—n,n + 24, ¢,| 
Tle) 


for large n. For an extensive treatment of the asymptotic properties of these poly- 
nomials for large n, see [9]. In general, the convergence properties of (2.3) are 
superior to the original series definition if the polynomials 


—n,n + 27, ¢,- 
F, ( ’ ? w) 
+2: +1 3 + y d, 


are interpolatory or nearly so. For example, if r = s, w = 1 and ¢; — ¢; is never 
an integer or zero, for i # j, 


—n,n + 2y, ¢, ~~ = 2c, |. 1 
to) [ie eae) fi 0()] 


p=DLa-DL ad, 


t=1 t=1 


(2.15) 


and the A and B,’s are constants independent of n. Incorporating the inequality, 
see (1), page 49 of [10], 


(2.16) \J(z)| s 





— exp {Im(z)}, 


we see that (2.3) converges like 


(2.17) 2 ta +2 a} lz". 


n=l 
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The original series definition, however, converges like 


(2.18) > ¥ zI". 


3. Further Expansions Involving Free Parameters. If in (1.3), we replace 
p by p — 2 and set a = a,4, 8 = ap, w = O, we get 
— Ss (ay) n(—z)” ( n+ a, | ) 
(3.1) 1 = 2 ®t paeni et n+b,,2n+7+1/\" > 


Then replacing a, by a, + k, b, by b, + k, y by y + 2k and multiplying both sides 
of (3.1) by 





(¢,)e(2w)* 
(d.)x k! , 





we get 


(cr)x(zw)* _ (¢,)e(2w)* = (k + ap)n(—z)" 








(d)ek! = (d,)e kl =a (kh + by)a(n + 2k + y)an! 
n+k+a, R 
x oFea(, k+by,2n+2k+y+1 :) 
(3.2) 





_—S& (ay)n(—2z)" ( a+6, | ) 
~ 2. Gar + ment Hn + b,,2n ty 44/7 


ye (= mala + v)e(ba)e(cr)e(w) 

(ap)e(d,), k! 

Summing the terms represented in (3.2) from k equal zero to infinity, and inter- 
changing the summation processes with respect to k and n, we obtain 


Cr | al “ (ap) n(—z)” ( n + ap | ) 
a (;;|=«) = x (ba)n(¥ + 2)an! Met \y + b,,2n+7+1|" 


—n,n T= 
x atrtel pie ( et . w) ° 
Pp» “s | 








(3.3) 


Using the Laplace transform techniques of Section 1, we arrive at the expansion 
Cr, et | s+ (et) n(ap)n( —2)” 
a ; aw) = ud 
«ea | d,, ful x (fu)n(ba)aly + )n”! 


Nn + dy,n + & | 
ted X arlerets (, + ban + fun ty + il) 


—n, nN 6 9 Cr 
X etrsaF oe ( T ¥5% wv). 


a,, ds 





Equation (3.4) and its confluent form not containing y are generalizations of results 
given by Chaundy, see (11), page 187 of [2]. 

As a final example of how Laplace transform techniques may be used in general 
expansions, we prove 
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sa(fe|r0) at (Gj | ne) = QO" 


q n—o (6,),n! 


pos — _ fe p+eti 
X ose ( n, 1 ? mer fe| ie), 


(3.5) 


1—n-—a,,d, | NN 
Again the proof proceeds by induction on p, g, r and s. The case p = gq = r= 8s = 0 
reduces to 
ete - . (r + hu) z 


n=0 n! 





(3.6) = (yz)" 
a ee 
<< n! Fo( a “), 
since 
(3.7) WF o(a|z) = (1 — z)™. 


Assuming (3.5) true for p, g, r and s, the inductions with respect to r and s can 
be effected by using with respect to » the Laplace transform techniques illustrated 
in the proof of (1.3). The inductions with respect to p and g are similarly carried 
out with respect to \ after making use of the relationship 


—n,@|_\ _ (a)n(—z)” —n,1 — n — b,|(—1)*” 
G8) wake ( by ?)  — uP ( l—n-4a, fo". 


This completes the induction proof. Equation (3.5) is a generalization of results 
given by Chaundy, see (12)—(15), page 187 [2]*. 
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TECHNICAL NOTES AND SHORT PAPERS 


Evaluation of Artin’s Constant and the 
Twin-Prime Constant 


By John W. Wrench, Jr. 


1. Introduction. In 1914 A. J. C. Cunningham [1] investigated the number of 
primes having the integers between 2 and 12 (exclusive of powers) as primitive 
roots. In particular, for the integers 2 and 10 he tabulated these counts for each 
of the first ten myriads, and deduced empirically that the density of such primes in 
the set of all primes in a given interval lies, effectively, between 0.355 and 0.400. 

According to H. Bilharz [2], E. Artin in an oral communication to H. Hasse on 
13 September 1927 proposed the question, whether, corresponding to a given non- 
zero integer a, the set of all primes of which a is a primitive root possesses a density. 
Artin’s conjecture as stated by Hasse [3], namely, that corresponding to every 
integer not a square and different from —1, there exists an infinitude of primes 
having that integer as a primitive root, has never been proved. Furthermore, Artin 
conjectured that the density of such an infinite set, if it exists, is the same for all 
non-power integers, and is given by the convergent infinite product 


x ee 
4 m1 p(p — mE 


taken over all primes p. This we have termed Artin’s constant. 
The twin-prime constant is the first member, C, , of a set of constants defined in 
1923 by Hardy and Littlewood [4] by means of the relation 


c= O44) G=3}- 


In that paper (p. 44) Hardy and Littlewood conjectured that 





" dz 
Pa(n) ~ 20+ [ (log x)?’ 


where P,(n) is the number of prime-pairs less than n, and they gave values of the 
two sides of this asymptotic relation when n = 10°(10°)10°. 

Apparently unaware of this paper, Sutton [5] in 1937 studied the average distri- 
bution of twin primes, using a probabilistic approach, and presented detailed em- 
pirical evidence (based on counts to 8-10°) for the validity of the foregoing asymp- 
totic formula. Subsequently, C. R. Sexton [6] performed an independent count of 
prime-pairs less than 10°, and revealed discrepancies in Sutton’s data as well as in 
the counts made by several other workers in this field. 

The most extensive and reliable empirical knowledge of this kind that is avail- 
able at the present time appears in an unpublished table of D. H. Lehmer [7]. 
Included in this table are the cumulative totals of prime-pairs in successive millions 
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as far as 37-10°, as well as the corresponding values of 2C2Ji.(x), where li,(x) 

represents [ (In t)~? dt. For example, Lehmer gives 183,728 as the count of prime- 
2 

pairs less than 37-10°, and gives 183,582 as the corresponding value of 2C2li,(z). 


2. Evaluation of Artin’s Constant. The present calculation of A has been based 
on the equivalent representation: 


In A = ¥,in(1-2) + yin(i-8 a En(i-2), 
p22 Pp p22 P p22 P 

where a and # represent the zeros of the polynomial p’ — p — 1. Hence, if the loga- 

rithms are expanded in Maclaurin series and if Newton’s formulas are used to 

express the sums of the same powers of a and £6 in terms of the coefficients of the 

polynomial p’ — p — 1, there results the series 


—nA = 34 p’+iud p +ia dp +---, 
p22 p22 p22 
where a, = a); + a2 + 1, fork = 2, and a = 0, a, = 2. 


The calculation was greatly expedited by computing separately the product of 
the first eleven factors of the product defining A, and then applying the preceding 
transformation to the infinite product consisting of the remaining factors. 

Carefully prepared tables of ) p * to 50D for k = 2(1)167 have been pub- 
lished by R. Liénard [8]. Subtraction from these data of independently computed 
and checked sums of reciprocal powers of the first eleven primes to 55D yielded 
values of a p * to 50D for k = 2(1)28. Direct summation was used to obtain 
these sums corresponding to k = 29(1)35, to at least 54D. Finally, multiplication 
by the appropriate coefficients a,/(k + 1) yielded approximations to the individual 
terms of the modified series for In A that were correct to at least 46D. 

Accordingly, the following rounded approximation to Artin’s constant is believed 
to be correct to 45D: 


A = 0.37395 58136 19202 28805 47280 54346 41641 51116 29249. 
3. Evaluation of C,. The present calculation of the twin-prime constant, C2, 


was performed in a manner similar to that of Artin’s constant. 
Since the product representation of C, can be written in the form 


2f g Cont Bel 
c= T{1- G45}. 


we have 
—In C, = 4, Dp? + 4h Dp +h Dp t+-::-, 
p>2 p>2 p>2 


where b, = 2*** — 2. 

In this calculation the product of the first ten factors was first found directly, 
and then the previously computed values of Deen? were combined with the 
appropriate coefficients b; to yield the terms of the modified series. 
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The resulting value of C, truncated at 42D is: 
C2 = 0.66016 18158 46869 57392 78121 10014 55577 84326 23 ---. 


This approximation confirms the accuracy of the 10D values given by Lehmer 
[7] and by J. C. P. Miller [9]. Sutton’s value of 1.3202 for Cy (= 2C2) is too low by 


about a unit in the last place. Recently Fréberg [10] has determined C, to 10D, of 
which the first 8 decimals are correct. 


As a by-product of this calculation we deduce a new approximation to D.,, 
defined by Rosser [11] as lim,.... D, where D, is determined by the relation 


1(-2)- <2 
Il 1-—-—j]= , 
m=? Pm (log n)? 
where p,, represents the mth prime. 


The value of D.. can be deduced from that of C, by virtue of the relation [9] 
| = 4e™*? C, : 





where y is Euler’s constant. 


Several years ago the writer computed unpublished values of both e” and e ” to 
170S. Consequently, the approximation of e°’ to about 50D was easily accom- 
plished, and the resulting value of D.. to 40D was found to be: 


De = 0.83242 90656 61945 27803 08059 43531 46557 50462 - -- 


This confirms the accuracy of the 12D approximation found by Rosser [11]. 

The writer should like to acknowledge the assistance of Dr. Daniel Shanks in 
the search for relevant literature and for several constructive suggestions in the 
preparation of this paper. 
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Iterated Square Root Expansions for the Inverse 
Cosine and Inverse Hyperbolic Cosine 


By Henry C. Thacher, Jr. 


Abstract. Let Ri = 1/2 + 27, Rin = 724+ R,. Then2 /(2—R,| and 
2'{|6 — 2+/3 + 3R, |}'” both converge to arccos z if | z| < 1 and to arccosh z 
if 1 < x < ~. Truncation errors for the two expressions are of the order of 2.” and 
2 respectively. 


1. Introduction. The availability on several modern automatic computers of 
square root operations which are approximately as fast as multiplication or division 
encourages investigation as to the uses which may be made of this operation in 
computation. Hammer [1] has described an iterative procedure based on square 
roots for finding cube and other odd roots which converge more rapidly than the 
customary iterations, but very few other authors appear to have considered this 
problem. It is the purpose of this contribution to describe a set of rapidly converging 
square root expansions for the inverse cosine, inverse hyperbolic cosine, and hence 
for the natural logarithm. 


2. Derivation. We start from the familiar identity 
(1) o(xz) = 26(3 V/2 + 2z) 


where ¢(z) denotes either arccos x(—1 S x S 1) or arecosh z(1 S zg S @). We 
restrict the multiple-valued inverse cosine to the branch (0 S arceos z S zn), 
the inverse hyperbolic cosine to the positive branch (0 S arccosh z), and take all 
square roots positive. Applying (1) repeatedly, we have 








(2) (x) = 2(26(3{2 + V2 + 2z}™)) = 2o(4{2 + V2 + 22})” 
(3) = 29(3{2 + [2 + V2 + 22]'7}"”) 


and, after k applications, 


(4) o(x) = Mo(a{2 + [2+ --- + V2 + 2a}'7}™) = o(4Rs) 


where R, contains k square roots. 
We may observe that R;, is an increasing function of z, and that for z = 1, the 
innermost square root becomes equal to +/4, so that for any k, 


(5) /2 < R(x) <2 (-l1 <2 1) 
and 
(6) 2 = R,(2) (ls2z< oo). 


Furthermore, since for —1 S z S 1, +/2 + 2x = 2z while the contrary is true 
for xz > 1, R,(z) is an increasing function of k for | x | < 1, and a decreasing func- 
tion of k for z > 1, and thus approaches 2 as k increases. Although we have only 
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proved this limit for real z > —1, it can be shown to be true for all finite real or 
complex z. 


If we multiply both sides of (4) by 2“, and take the cosine or hyperbolic cosine, 
we find 


(7) ¢ (2 “o(x)) = 4Re. 


Now both the cosine and the hyperbolic cosine may be expanded in Taylor’s 


series. If we use only the first two terms of the series, and the error term, EF; , we 
find 


s) 1 a (0) 


where the upper sign is to be used for the inverse cosine. If we use the first three 
terms, we have 


+E, =5R 


—k 4 
(2 ¥)) + E; =-1R,. 


Solving (8) for ¢(x), remembering (5), (6), and our restrictions of the values of ¢ 
to positive quantities, we find 


a 





() 1 4 4) 








(10) o(z) = 2 +//2— R, + 2B; |. 
Equation (9), a quadratic in (2“¢(2) )’, has the root 
(11) (2"o(z))’ = |6 — 2/3 + 3R, — 6B, |. 


The second root is easily seen to give a value for ¢ outside the specified range. 
(For |x| < 1,2"(6 + 2+/3 + 3R, — 6E;) > 2"(6 + 2{3 + 3 +/2) }"" > 22’, 
while for zx < 1, —6 — 24/3 + 3R:. — 6E; < 0). Hence, 

(12) o(x) = 2{|6 — 2/3 + 3R. — 6E;|}™. 


The desired approximations are, of course, (10) and (12) with the error terms EF; 
and E; omitted. 











3. Truncation Error. Inverse Cosine. In estimating the truncation error incurred 
by neglecting £; or EZ; , it is more convenient to analyze the expansions for arccos x 
and arccosh x separately. For | x | < 1, (10) becomes 


(13) arccos x = 2°4/2 — R, + 2B; = 40/2 —Ri tn 

while (12) becomes 

(14) arccos x = 2°{6 —2 4/3 + 3R, — 6B;}*” = 2{6-2 V3 4+ BR. } 4+ ws. 
By Taylor’s theorem, expanding (13) in powers of E; , 

(15) arccos r = 2*{+/2 — R, + 1/+/2 — R, + 208; E;} © < 6 < 1). 


Now the error ir the cosine expansion (8) is of the same sign as, and less in magni- 
tude than the first neglected term, so that 


—k 4 
(16) 0<h< (2 eens #) 





1e, 


bd 


rs 
we 


ee 


5 
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Hence, 

(17) 1/2~/2 — R, = 1/2~V/2 — Ri + 208; 
and 

(18) 0 S m < 2“ (arecos z)*/24./2 — Ri 


and the upper bound on 7; is of the order of 2 (arccos z)*. Expanding (14) in the 
same way, we find 


arcos z = 2° {6 — 23 + 3h} 
(19) + 3E;/[V3 + 3R. — 66Es {6 — 2/3 + 3k. — GOEs}""| (0<6@<1) 
= 24{6 ~ 20/3 + 3h} + 





while 
: (2~ arccos x)° 
> ME 
(20) O2>E=2 720 
so that 
—4k 5 
(21) (ze2 —* (arecos x) 


ee S200 9/3 + BR 
and 7; is of the order of 2™. 


Thus, our two approximations have errors of opposite sign, and provide bounds 
on the true value. 


4. Truncation Error. Inverse Hyperbolic Cosine. For the inverse hyperbolic 
cosine, x 2 1, and we have: 
(22) arecosh x = In (2 + +2? — 1) = 2/R, — 2 — 26, = 2VR—-24+ m 
(23) arecosh z = 2*{24/3 + 3R. — 6E; —6}'" = 2° {24/3 + 3R, —6}"" + ws. 
Again using Taylor’s theorem, 
(24) o/h, 2 = = 2 Vi = 2 - - ge} (0<0<1). 


The remainder in the hyperbolic cosine series is 





(2~“areccosh x)* 











(25) E; = cosh (2~*6’ arccosh x) a (0 < @ <1) 
and has bounds 

(26) O<sE<a cosh a ) (2~* arccosh x)‘ 

so that 

ry > qe 2‘E; . 2*E; 2‘E; 


JR, -2—-08, VPs —-2—Es ~ 2-* arecosh z* 
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Hence, 7; has the bounds, 





(28) 02m z= — 2 errcomh) 

















ie _ 5% R, (arecosh x)* 
7A cosh (2° arecosh +) = —2 oa - 
In similar fashion 
3Es 
(29) m= = —2" seis ites 
: (2 /3 + 3R, — 6B, — 6}"” +/3 + 3R, — 6OEs 
> — 2 3Es 
(30) {2-V3 + 3R. — 6E; —6}"" V3 + 3R: — 6Es 
—— 6E; 





(2-* arecosh x) [6 + (2-* arccosh x)*]" 
The remainder £; is given by 

















g-*, L ~)\6 
(31) Es = cosh (2~*@’ arccosh x) (2 some 5) (0< @ <1) 
so that 
—k J 

(32) O<E< cosh (2 — ) (2~ arecosh x)° 
and 

ie os oe cosh (2~ arecosh zx) __2™(arecosh az) be 
(33) AA 120 [6 + (2- arecosh x)?] 


t: 7° &. (arecosh x)° 
6 + (2— arecosh x)? 240 1 





The error bounds given do not appear to be unduly conservative, and are ap- 
proached by the actual error as k increases. 


5. Roundoff Errors. These algorithms are subject to serious roundoff error when 
k is large, and, except for special cases (e.g. z = 1) are incapable of evaluating the 
functions to the full accuracy of the arithmetic being used. This is so with either 
fixed or floating arithmetic, since for both algorithms, 2~“¢'(x), a relatively small 
quantity, must be calculated from the difference of two quantities which are each 
greater than one. Actual computing trials have indicated, however, that approxi- 
mately three-quarters of the total number of digits carried may be obtained cor- 
rectly, even when particular attention is not devoted to minimizing roundoff. 


6. Discussion. The algorithms have three major advantages: (a) they are simple 
to program (and rapid to calculate when automatic square root operations are 
available); (b) they require only one or three stored constants; and (c) they have an 
extremely wide range of acceptable convergence. As can be seen from the error 
bounds, these expansions converge over the entire range —1< x < ©, and the con- 
vergence especially for (14) and (23) is notably rapid compared to the power series. 
For example, using a programmed 16-significant-decimal digit double precision 
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interpretive routine for the LGP-30, and (23), arccosh 250.001 (i.e. In 500) was 
found correct to 4 decimal places with k = 5 (6 square roots) and to 11 places with 
k = 10 (11 square roots). 

The convergence of the sequence of approximations is only first order. At some 
cost in programming effort, it would be clearly possible to increase the convergence 
by one of the standard extrapolation techniques for accelerating the approach to 
the limit. However, the excellent convergence already present makes it unlikely 
that this device would be worthwhile unless the need for high accuracy was such 
that it was essential to keep k as low as possible. 

A special case of (13), with 2 = —1, has been known for a long time. This 
expansion, which has the limit +, can be obtained as one-half the perimeter of a 
2*-gon inscribed in a circle of unit radius [2]. However, the general case, and the 
expansions obtainable by retaining the fourth-degree terms in the series for cos zx 
and cosh x appear to be new. 


Argonne National Laboratory 
Argonne, Illinois 


1. P. C. Hammer, ‘Iterative procedures for taking roots based on square roots,’”’ MTAC, 
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An Eigenvalue Problem Arising In Mass And 
Heat Transfer Studies 


By J. S. Dranoff 


1. Introduction. In a recent paper [1], 8. Katz has considered the problem of 
catalytic chemical reactions occurring on the inside surface of a cylindrical tube. 
For the case of laminar flow of reactant through such a tube, he has shown how one 
may generate basic kinetic data for the reaction in question from easily made over- 
all conversion measurements. The interested reader is referred to the original paper 
for the details of this analysis and its application. 

In order to make use of Katz’s analysis, one must have on hand the solution to 
the following Sturm-Liouville type eigenvalue problem: 


a (= r) 


- = 
dx\" dx ) +r 4e(1 © )on(x) = 0, 0 


lA 
& 
IIA 


(1) o.(x) regularat zx =0 


on'(1) = 0 


where the ¢,(x) are the eigensolutions and the X, are the eigenvalues, with n = 
0, 1, 2, --- . The first boundary condition leads, as in the case of Bessel’s functions, 
to the condition ¢,’(0) = 9. 
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In addition to ¢,(z) and X, , it is also necessary to have available the numerica: 
values of the following integrals in order to complete the analysis. 





(2) a, = | & (a axl — 24) de 
N, 
(3) p, = [[ 220(2) ar 
N, 
(4) Ni = [ on (x)4x(1 — 2°) dz. 


The reader may recognize that a, and b, are the coefficients in the expansion of 
the functions 2’/2 and 2z/4x(1 — 2’), respectively, in a series of the eigenfunctions 
¢.(), with a weighting function of 42(1 — x’). Further, the integral N,, may be 
regarded as a normalization factor. 

It should be noted that problems similar to (1), but with alternate boundary 
conditions, have been studied by many investigators in the past. The earliest work, 
which dates back to 1883, seems to be that of Graetz [2], who considered the case 
in which the second boundary condition of (1) is ¢,(1) = 1. Most recently, Brown 
[3] recalculated the solutions to Graetz’s problem by numerical techniques using a 
digital computer. However, no solutions pertinent to the present boundary condi- 
tions are available. Hence, the solution to the problem (1) and the determination 
of the integrals (2), (3), and (4) was the objective of this work. Since no complete 
analytic solution was possible, it was decided to solve the problem by numerical 
means using a Burroughs Datatron 205 digital computer. In addition, an effort 
was also made to obtain an asymptotic analytic representation of the solutions to 
the eigenvalue problem. 

The normalization ¢,(0) = 1.0, for every n, was adopted for convenience in 
numerical calculations, rather than N, = 1, as chosen in [1]. Examination of the 
defining differential equation in (1) shows that any eigensolution may be multiplied 
by an arbitrary scale factor and still remain a solution. Thus, it is perfectly legiti- 
mate to make this normalization. 


2. The Asymptotic Solution. An asymptotic solution to the problem (1) with 
the specification that ¢,(0) = 1.0 was attempted for large n. The method used was 
exactly that used by Sellars, Tribus and Klein [4] in dealing with the Graetz prob- 
lem. Their procedure was followed except as required by the different boundary 
conditions of (1). The result was a three-piece analytic approximation to the eigen- 
function ¢,(x) and an equation for estimating the eigenvalues \, . 

The actual solutions obtained were as follows: 





(5) n(x) = Jo(2a~/rn) for x near 0 
¢.(z) = 1 cos [x(1 — 27)"?4/X, + Vd, aresin x — 7/4] 
(6) 4 WAL 24 — 2)" 


for 0<2< 10 
(7) gn(z) = (—1)"*V3(1 — 2)? J a(3(1 — 2)*?V/82,) for x near 1.0 





f 
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where the eigenvalues may be approximated by the expression 
(8) An = 4(n + 1/3)’. 


These eq. \tions are valuable for checking the results of subsequent numerical 
solutions. In particular, (8) is most useful in predicting the approximate location 
of the eigenvalues \,, . 

Although it is difficult to predict the accuracy of these asymptotic expressions, 
one may expect them to be quite good for n = 5, judging from the results of Sellars, 
et al. 


3. The Numerical Solution. The numerical solution of (1) was carried out by a 
trial-and-error procedure. A value of \,, was assumed, and the differential equation 
was integrated numerically from z = 0.0 to z = 1.0, using a standard Runge-Kutta 
method [5]. The initial conditions ¢,(0) = 1.0 and ¢,’(0) = O were assumed, as 
discussed earlier. If the resultant solution satisfied the criterion that ¢,’(1) = 0.0, 
it was deemed an eigensolution of the problem and the corresponding A, a true 
eigenvalue. If this condition was not satisfied, a new value of \, was chosen and the 
process repeated until successive values of \, agreed within six (6) significant 
figures. (It should be noted that all computations were done in automatic floating 
point form using ten (10) decimal digit words, consisting of a two (2) digit charac- 
teristic exponent and an eight (8) digit mantissa. ) 

Successive estimates of \, were made according to the following formula: 


(9) rt = a, — 2's 


where the superscript indicates the “ith” iteration and the value of m is given by 


(10) ~~ [¢,.’(1)] a [o.’(1)) 


@ =D 
An — An” 





Equations (9) and (10) represent a simple linear approximation to the ¢,’(z) — Aa 
relationship near z = 1.0. Since this method requires the results of two previous 
iterations in order to estimate a new value of X, , it could not be applied until the 
third iteration. The second approximation to \,, was therefore found by arbitrarily 
increasing the starting value by 0.5. By using the predictions of equation (8) to 
estimate \,”, it was normally possible to make the calculations converge within 
four or five iterations. 

A fixed interval size was used in each integration. The final value was chosen 
such that a smaller interval would not produce any change within six (6) significant 
figures in the calculated eigenfunctions and eigenvalues. The number of intervals 
was varied in the calculation of successive eigenfunctions, but was always greater 
than 25 n. This criterion was established in order to provide for the increasing 
complexity of the ¢,(z) within the interval [0, 1] with higher values of n. 

When the eigenvalues and corresponding eigensolutions were obtained, the 
integrals (2), (3), and (4) were evaluated by numerical integration on the com- 
puter, using a standard Simpson’s rule integration procedure. 


4. Results. The numerical solution described above was carried out for n from 
1 to 20. Table 1 presents the results of these calculations in terms of the eigenvalues, 
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TABLE 1 
Numerically Calculated Eigenvalues, Eigenfunctions, and Derivatives 





An @» (1) on’ (1) Ax* 














n | 

1 6.41990 ~0.492517 | —1.0 x10-" | 0.02 
2 20 .9654 +0.395509 | -4.0x10-° | 0.0100 
3 43.5417 —0.345874 | +2.3xX 107 | 0.0050 
4 74.1341 +0.314047 | —3.7 x 107 0.0050 
5 112.737 | —0.291253 | +1.5 x 10-6 0.0050 
6 159.347 +0.273810 | ~—1.5 x 10-° 0.0050 
7 213.963 —0.259852 | +3.5x10-> | 0.0025 
8 276 .583 +0.248329 | -7.5X10* | 0.0025 
9 347 .206 —0.238591 | +4.5 x 107 0.0025 
10 425.835 | +0.230189 | —2.5 x 10+ 0.0025 
11 512.462 | 0.222865 | +5.0 x 107 0.0025 
12 607 .093 | +0.216373 | -9.6x 10% | 0.0025 
13 709.725 —0.210568 | —7.6 x 107 0.0020 
14 820 . 360 +0. 205334 4+5.3x107 | 0.0020 
15 938 .997 —0.200580 +1.2x107 | 0.0020 
16 1065.63 +0. 196234 —5.5 X 10-* 0.0020 
17 1200.27 —0.192237 | +42.1x10-> | 0.0020 
18 1342.92 +0.188546 | -3.2x10-* | 0.0020 
19 1493 .56 —0.185120 +1.5X10-° | 0.0020 
20 1652.20 | +0.181929 | +1.8x10-* | 0.0020 





* The number of steps used in the numerical integrations is equal to 1.0 Ar. 











TABLE 2 
The Coefficient and Normalization Integrals 

n | an bn | Nn 

1 —0.243568 — 1.02664 +0.190138 
2 +0. 124300 +1.06525 | +0.107731 
3 —0.080695 —1.08411 | +0.075286 
4 +0.0585046 +1.09553 | +0 .0578066 
5 —0.0452500 —1.10329 | +0.0469413 
6 +0.0365248 +1.10895 +0.0395161 

7 —0.0303940 —1.11331 +0.0341194 
s +0.0258748 +1.11677 +0 .0300203 
9 —0.0224156 —1.11957 +0.0268008 
10 +0.0197100 +1.12198 +0 .0242052 
11 —0.0175137 —1.12392 +0.0220681 
12 +0.0157222 +1.12564 +0 .0202777 
13 —0.0142302 —1.12714 +0.0187560 
14 +0.0129707 +1.12847 +0.0174469 
15 —0.0118950 —1.12964 +0.016309 
16 +0.0109669 +1.13069 +0.0153098 
17 —0.0101598 —1.13165 +0.0144264 
18 +0 .00945062 +1.13250 +0.0136393 
19 0 .00882477 —1.13229 +0.0129337 








20 +0 .00826766 +1.13400 +0 .0122975 
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and the eigenfunctions and their derivatives evaluated at = 1.0. In addition, the 
size of the integration interval used in the numerical calculations has also been 
tabulated for reference. Note that the leading eigensolution for n = 0, namely 
do = 0, do(x) = 1, has not been tabulated. 

More extensive values of the eigenfunctions are shown in Table 4 for intervals 
of x = 0.1. Complete tables are available through the Engineering Research Section 
of the American Cyanamid Company, Stamford, Connecticut. 

The integrals (2), (3), and (4) were evaluated for each of the above functions 
and are presented in Table 2. The data of Tables 1 and 2 may now be used to 
evaluate the function 


(11) M(0) =1+ > Onna + Ba)da(1e™* 


which is equation (48) of [1]. 

The values of a , by , and No have not been tabulated. These may be readily 
seen to equal 3, 1, and 1, respectively, by examination cf equations (2), (3), and (4). 

Some idea of the validity of the asymptotic expansion described above can be 
obtained by comparing the calculated and predicted vaiues of the eigenvalues as 
determined by numerical integration and equation (8), respectively. These data 
are shown in Table 3. It is apparent that the asymptotic expression for the eigen- 
values is quite a good approximation for n = 5. It is assumed that this agreement 
holds for the asymptotic forms of the eigensolution as well, although no specific 
test has been made of this point as yet. 














TABLE 3 
Test of the Asymptotic Expression for the Eigenvalues 
”" d (Calculated) A (Predicted by Equation 8) Apred/Acal 
S. = 
1 6.41990 | 7.11111 0.9028 
2 20 . 9654 21.7778 0.9627 
3 43.5417 44.4444 0.9797 
4 74.1341 | 75.1111 0.9870 
5 112.737 113.7778 0.9909 
6 159.347 160.444 0.9932 
7 213.963 215.111 0.9947 
8 276.583 277.77: 0.9957 
9 347 .206 348.444 0.9964 
10 425.835 427.111 0.9970 
11 512.462 513.778 0.9974 
12 607 .093 608 .444 0.9978 
13 709 . 725 711.111 0.9981 
14 820.360 821.778 0.9983 
15 938 .997 940 .444 0.9985 
16 1065 .63 1067 .11 0.9986 
17 1200.25 1201.78 0.9987 
18 1342.92 1344.44 0.9989 
19 1493 .56 1495.11 0.9990 
20 1652.20 1653 .78 0.9990 
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The recent work of Siegel, Sparrow and Hallman [6] has just come to the 
attention of the author. These workers have considered this problem in the heat 
transfer context. They report values of the eigenfunctions and eigenvalues which 
are in excellent agreement with the more extensive data of the present work. 
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Efficient Continued Fraction Approximations 
To Elementary Functions 


By Kurt Spielberg 


1. Introduction. This paper describes an application and extension of the work 
of H. J. Maehly [1] on the rational approximation of are tan z, and of E. G. Kog- 
betliantz [2], who developed Maehly’s procedure so as to be applicable to the com- 
puter programming of elementary transcendental functions. 

It is to be shown here that certain modifications, such as the introduction of 
terms which are easily computed on specific computers, lead to considerable im- 
provements. In particular, the application of the modified method to several ele- 
mentary functions will be described and corresponding final results will be given. 
Some of these approximations have been used with great success to develop sub- 
routines for the IBM 704 and 709 computers. Our experience indicates that the 
method of Maehly and Kogbetliantz, as modified below, is superior to other current 
numerical procedures. 


2. The Modified Method of Maehly and Kogbetliantz. The basic idea made use 
of by H. J. Maehly in connection with f(z) = arc tan z is to approximate the func- 
tion f(x) by a ratio of two Chebyshev sums of order k 
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f(a) = Ya,-7(2) / [1 +dd. r.(2) 
q) 


+ az) / [1 + 2d b- 12) |= $*(2) +A 


where H(z) = >-%0A;”-Tx414;, and A is the absolute error. 

If the coefficients b, are small, which will normally be the case for —1 S x S 1 and 
reasonably rapid convergence of the power series for f(z), then the denominator of 
the error term is close to one over the interval of approximation and H(z) repre- 
sents the absolute error A with sufficient accuracy (compare [3]). The order k is 
chosen so as to keep A below the desired upper limit of accuracy. In order to 
evaluate the unknown coefficients a, and b, , one must know the coefficients c, of the 
Chebyshev expansion of f(x). A comparison of that expansion with (1) leads, after 
use of the identity 


(2) 2T n(x): T,(2x) = T min(Z) 7 T m—n(2) 


to a set of 2k + 1 simultaneous linear equations in the 2k + 1 coefficients a, and 
b, . Additional equations can be established for the coefficients A;“) in the error 
function H(z). For even and odd functions the subscripts of (1) are changed as 
follows: even, r — 2r, s — 2s, 7 — 2j, k — 2k + 4; odd, r— 2r + 1, s > 2s, j — 2), 
k— 2k + 1. 

When this scheme was applied in practice, several additiona! ideas suggested 
themselves. They can be listed briefly as follows: 

a) Application of the method to functions that can be expressed as ratios of 
Chebyshev series, such as tan (}72). 

b) Use of different degree numerator and denominator polynomials in (1). 

c) Consideration of unequal intervals for two complementary expansions, such 
as sin ax and cos Bz, a + B = 2/2. 

d) Reduction of the relative error by means of a linear correction term in a 
neighborhood of x = 0. 

e) Reduction of the error term through introduction of a new parameter that 
does not lead to a full additional multiplication. 

The first three points should become clear in the sequel and need little amplifica- 
tion. Point d is usually of concern for odd functions f(x), if it is desired to obtain 
accurate results for g(x) = f(x)/x as x — 0. Chebyshev methods applied to the 
function f(z) produce an approximation f*(z) such that the absolute error 
| f(x) — f*(zx) | is (approximately) minimized over an interval such as —1 S z S 1. 
The relative error of such an approximation, | f* — f |/f, usually becomes intolerably 
large as x and f approach zero. The natural way to cope with this difficulty would 
be to apply the Chebyshev approximation method to g(x) rather than f(x). Then 
the relative error |f* — f |/f = | x-g* — x-g |/x-g = | g* — g |/gis nearly minimized 
over the interval if g(x) does not vary too much. This approach, however, has the 
drawback that the improvement in the neighborhood of zero is paid for with a 
decrease in accuracy in the remainder of the interval. We have found that computer 
subroutines can easily be written so as to use f*(z) in most of the interval and 
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f*(x) + x-C in an appropriately chosen neighborhood of zero. The choice of the 
correction term C' will be discussed further below. 

The most important modification of the method, point e, arises if the special 
machine characteristics of digital computers are taken into account. The reduction 
of the error clearly depends on the introduction of the parameters a; and b; . Each 
new parameter allows reduction of one more error term A; to zero, but also results 
in an increase of the number of multiplications (or divisions), M, by one. We can, 
however, achieve a compromise by restricting the last parameter to a set of values 
which may allow the correspondingly introduced multiplication to be performed in 
a manner particularly suited to the calculator in question. In the case of the IBM 
704 the multiplication might be reduced to a shift (a “cheap” multiplication with a 
power of two), in the case of the IBM 709 to a variable length multiplication requiring 
little time. Among the permissible values for the newly introduced parameter, that 
one is chosen which allows a maximum reduction of the dominant error term. In 
other words, one of the residuals in the system of linear equations for the a; and b; 
is not reduced to zero but only below a certain value determined by the desired 
accuracy. 

As an example, we may consider as our newly introduced parameter the coeffi- 
cient as in 


(3) f(x) = (a, T, + a; T; + as Ts)(1 + b. T:)* = .f Ki+K2- 2+ P me ) 
4 


24 K. 
Evidently K. = 8a;/b2 . In accordance with the above discussion we restrict K, to 
the form 2” (n --- any integer), so that a; becomes restricted to the set of values 


b.-2”. The integer n is chosen so as to reduce the absolute error as far as possible. 
Numerical details will be given in Section 3. 

The modified procedure of Maehly and Kogbetliantz can now be outlined for- 
mally as follows. Given a function f(x), find an approximation /*(z) 








20 I 
> «: T(x) ) Bs a; T;(x) 
(4) f(z) = =>———__, f(z) = ——= 
D 4; T(x) 1+ 0b: T(x) 


To determine the coefficients a; and b; , we consider the absolute error 


f(z) — f*(2) = [= aT; (1 +d, r)| 


i=1 


i=l i=0 1=0 i=0 


” [= d; T; (1 + > bi r.)| : > R;- T(z) 


i=0 


(5) [+ den) Lean Lar Lar| 


The coefficients R; can be viewed as the residuals, usually rapidly diminishing in 
magnitude as 2 increases, of an infinite number of linear equations in the s unknowns 
ti: b:,be--+ bm, ao,a°** a. 
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(6) Ri= Lew asthe, i= 0,1,2,3,---,8-1, --- 
j=l 
The e;; and f; are simple sums formed with the coefficients of the given function 


f(z), c; and d; . For instance, the important special case of a polynomial f(x) gives 
rise to the following residuals: 


Ro = coo +4 Di bata — Go, by =1 
n=1 


(7) lsitsm: R =e + $0); + +> bn: (Cn4i + Cini) — a; 
1 


a= 


m <i: Ry = 4D0 ba-(engi + Cin) — G; 
n=0 
For iz > I, the term —a; is omitted above. 


Usually one sets the first s residuals equal to zero so that the final error is deter- 
mined by the absolute sum of the remaining residuals. Instead of this we endeavor 


to reduce the first s + 1 residuals below a desired bound 4, by introducing the 
additional parameter a;4, 


R; = w;-4, z= 0, 1, ---,s—l,8 
(8) 


Zou = Gin = k-2,, & yy 


The w; are weightfactors which can be chosen arbitrarily, usually as 0 for i < s 
and as 1 for i = s. The choice of z, and k depends on the transformation from the 
rational approximation to the continued fraction. In the example given above 
z, is equal to b. and k is chosen to be of the form 2”. 

The residuals R; clearly become linear combinations of s + 1 variables z; . 


In view of (8), however, they can be expressed in terms of the first s variables and 
k. 


lA 
IIA 


Vv 


(9) Re = wed = > e4;-4; + eseuk-t, + fi, 1=0,1,2,---,s—1, 1 8. 
j=l 


These equations can now be solved for the x; in terms of k. As a consequence, one 
can determine the residual R, as a function of k 


(10) zi=2(k), «= 1,2,---,8, R, = w,°6 = f(k). 


Finally one chooses among the manifold of permissible k, namely of those k which 
permit the replacement of a multiplication by a more favorable operation, that 
value which minimizes R, . It is usually possible to reduce R, so substantially that 
the leading term of the absolute error becomes R,4; . Except for a bounded factor 
stemming from the denominator in (5), the final absolute error is given by 


(11) AbD |wePs|+ Do |ReTi| © 2d | Ril. 


i=s+ i=s+ 


It is perhaps of interest to point out that, when applied numerically, this procedure 
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usually produced values of k which did not only reduce R, but also décreased R,,, 
in magnitude. 

We finally turn our attention again to the correction term discussed in point d. 
Inspection of (11) indicates that in a sufficiently small neighborhood of z = 0 one 
can approximate A by the leading term R,4:- 7.4; = R.4:-const-z. The size of the 
interval about zero will primarily depend on the relative magnitude of the two 
lowest degree terms in the replaced Chebyshev polynomial. For instance, 
| Ro- Ts| =| Re: (2562" — 576x” + 432z° — 1202° + 9z) | can be replaced by | 9-Re-x | 
for |x| << 3/+/120, or for a neighborhood of zero in which 120-Ry-z* < wy-6. (Of 
course, Ry-T, must also be below the tolerable error limit.) From a practical 
standpoint, it is simplest to compute the coefficient of the linear error term as 
lims+o [f(z)/z — f*(x)/z]. 


3. Applications and Results. The method outlined above was applied to produce 
rational approximations for the development of optimal elementary function sub- 
routines for the IBM 704 and 709 computers. The resulting routines were tested 
carefully and have been found to give, within the limits of roundoff error, results 
of the predicted accuracy. 


a) Sine Approximation, 8-Digit Accuracy. 


(12) 2+ Ky, 


sin* az = 2(a, 7; + a,7T; + asTs)(1 + b, T2)* =z (K+ 4a + we) 
a = 3, z= 32, —-352=53 

In this case (10) becomes 

R; = 10°°-[—.28029 + 10*-.38076(+2”"°-a*® + 10°°-.55934)~"). 


The minimum of R; is reached for n = —3. The correspondingly attained im- 
provement becomes apparent if one compares R;)_; with R;)_.. , the error without 
correction term. 


Rust SbxXwW”, MR. 4xw’ 


In the sequel we shall give the results of our computations as they were calculated, 
that is to more places than is usually warranted by the accuracy. 


a, = .14838 52081 2231 K, = —10° X .19845 92426 192 
a; = —10 ° .49269 95891 193 K;= 10° X .10429 26708 144 
a5= 10 °.37911 69631734 = 2”°-a'b, Ks= 10° X 50030 24548 541 
b = 10° .89864 76164 110 

a=.3, n=-—3, A (absolute error) > 10” X .55, R (relative error) 


& 14 A/.2955 = 10° X .26 
Check: 1) lim (sin* z)/z = .99999 999899 as z—-0 
2) sin* (.15) = .14943 81324 75 -- , sin (.15) = .14943 813--- 
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b) Cosine Approximation, 8-Digit Accuracy. 
cos* 1.32 = 2(aoT'> + a2T'2 + a47's + a67's)(1 + boT2 + O47 s)™ 
= K, a 22” + Ki{2 + K, + K;(z + Ks) J" 


z2 =132, —-13832513 

a = 30812 59625215 K,= 10° X .33947 71494 5237 
(13) a= —.17646 68549 891 K; = —10° X .29702 03659 0243 

a= 10° X 49379 28852647 K,= 10° X .57936 13928 3225 


a 
o 
Il 


—10* X .34496 21183611 Ks 10° X .14546 86265 7824 
10° X .20951 16328571 Ke= 10° X .48789 05740 6695 
bs = 10° X .81647 83866 535 
In this example equations (8) and (10) take the form: 
ae = 2”*-(1.3)?-b, Rs = 10 °(—.1490 + 2” .3608)(.2862 + 2” *1.667)* 
Rs)nex—o0 = —.52X 10°, Rs)nao& 20 X 10° 
A (absolute error) & .30 X 10°, R (relative error) ~ .20 X 10°. 


> 
I 


The actual computation of (13) involves the subtraction of two large numbers with 
a corresponding loss of accuracy. Hence a transformation to a more satisfactory 
continued fraction was performed 


cos* 1.32 = H, — 22 + (Hs + 3202*)[2 + Hs + Hs(2* + He) J" 

10° X .19477 14945 2366 Hs = 10° X .22874 43195 6870 
H; = —10‘ X .32763 39951 6402 H. = 10° X .24144 89469 4287 
H,= 10° X .82580 30199 5633 

Check: 1) cos* (0) = 1.00000 00005 
2) cos* (1) = .54030 23025, cos (1) = .54030 2306--- 


(14) 


c) Tangent Approximation, 8-Digit Accuracy. 
tanjrz = (a,T, + a7; + a57)(1 + beT')* 
(15) = 2[Ki + Ke + K3(2 + Kx) '] + (.165)10 z 
—_—_—_——” 


correction term 


—in S2=}07 Sh, correction term added for —.15 S z S .15 
a, = 86739 36410 Ki = .18717 82697 7 
a; = —10° .10096 94650 K:.= 10° .41503 90625 = 2’ (.100 01). 
a; = —10* 35876 64246 K;= —10 .20070 07228 1 
b= —  .1427391684 K,= —10 .24691 85502 1 


A (absolute error) & .842 x 10° 
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The relative error is reduced by the addition of the correction ternt. 
Check: 1) tan* z/z = 1.00000 00001 as z—0 
2) tan* (.1) -10033 4665, tan (.1) 

3) tan* (.7) .84228 83844, tan (.7) 


10033 467 --- 
84228 838 --- 


Il 
lI 


d) Cotangent Approximation, 8-Digit Accuracy. 
cot* inz = (aT; a a3T'3)(1 + beT 2 + bTs)* 
= 1/2[K, + Kee” + K;(2 + Ks)” — 526 X 10") 
Me 


correction term 


— 49 S z S ir, correction term added for —.15 S z S .15 
“ a, = 86369 96360 Ki= 10 X .34180 166678 
a; = —10" X .13359 274436 K:= — 10156 25 = —(.000 110 1). 
b, = .15019 244548 K;= 10° X .25226 53989 66 
b= 10° X 53281 462842 K, = —10° X .10432 74050 83 
A (absolute error) Y .263 X 10° to .86 xX 10° 
Check: 1) z-cot*z = 1.0000000525 — 526 x10" as 2-0. 
2) [eot*.1)’ = .10033 46678,  tan.1 = .10033 467... 
3) feot* 15° = .15113 5221, tan .15 = .15113 522... 


e) Tangent Approximation, 10-Digit Accuracy. 
tan* tnx = (a7; + a373 + a57's)(1 + T2 + biTs)* 
= 2{(Ki + Kez’) [2 + Ks + Ka(2 + Ks) ‘4 


Sz=i0trS}4 


ae 
qy 
IA 


a, = 86130 00805 00276 K, = —10' X .62993 45787 14378 
(17) a; = —10° X .15478 46841 31747 ‘K, = 06738 28125 
= (.10001 01).-2° 
a,= 10° X .23254 407227 K; = —10° X .17890 72313 8022 
b = — .15503 39460 54144 K,= 10° X .11511 65957 03706 


bh = 10° X .87881 25568 77435 K, = —10' X .99312 26653 90157 


It was again found desirable to transform to an equivalent approximation with 
smaller coefficients: 


tan* z = 2{(H, + Hee’) {2 + Hs + (Hs — 162°) (2 + Hs) T} 
H,=Ki, H:=K:, H;= —10' X .18907 23138 022 


(18) 


H, = 10° X .43783 03075 87186 H; = Ks 
A (absolute error) > .7 X 10", R (relative error) > .83 x 10°” 
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Check: 1) lim (tan*z)/z = 1.00000 00000197 as z—0. 
2) tan* }rz = (.78539 81635 2).4/z, ie = .78539 8163 -- 


f) Sine Approximation, 10-Digit Accuracy. 

sin* d4x = 2(a,T1 + a37'3 + as7’s)(1 + b2T2 + 47s) * 
2{K, + K.[z’ + Ks + Ki(z’ + Ks) J} 
2{H, + (Az + 62*)[2 + Hs + Haz’ + Hs) T"} 


—in Sz=}0t She 


a, = 36498 44708 84912 K,= 10' X .71483 02660 86945 
a; = —10° X .78595 99360 53327 K: = —10° X .80285 00024 48424 
a5 = 10° X .30425 2059202251 K;= 10° X .55072 65773 70549 

es b, = 10° X .10166 05194 0260 K,= 10° X .12558 99943 5548 
b, = 10° X .21677 076180 Ks= 10° X .16632 65254 62696 


H,= 10° X .11483 0266086945 H,= 10‘ X .21114 87369 0673 

H, = —10° X .37773 44209 59983 ss 

H;= 10° X .70852 59691 47400 

A (absolute error) 4356 X 10", _—&R (relative error)  .863 X 10°” 
Check: 1) lim sin* z/z = .99999 9999977 as z—-0 

2) sin* (.5) = 4794255386, — sin (.5) = .47942 5539 --- 


ll 


85271 33685 84500 


g) Cosine Approximation, 10-Digit Accuracy. 
cos* Sax = 2(aoT’o + a2T 2 + a47's + a6T's) (1 + b2T2 + yTs) 
= K, — 22+ Kifz + K;+ Ki(2 + Ks) J" 
= H, — 22 + (AH. + 3382")[2 + Hs + Hu(2 + Hs) Jy" 


—4in S2= 407 SS} 


a = 42553 53145 92886 K,= 10° X .33841 65629 20989 
a, = —10" X .69950 71904 10770 K. = —10° X .29473 89085 26762 
a, = 10° X .68476 48239 114832 K;= 10° X .57451 41742 44846 

(20) .ag = —10°° X .17046 17065 67264 K,= 10° X .14616 00034 97481 
b = 10° X .76660 47537 720 Ks = 10° X .48882 57695 11137 
bh = 10° X .11053 684400 
H, = .41656 29209 89 H,= 10° X .42283 05642 42266 
H,= 10° X .63340 59841 2301 H; = 39331 18492 483 


H;= 10° X .10594 06825 2635 
R (relative error) & 5 X 10™ 
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h) Logarithm Approximation, 8-Digit Accuracy. 
log.* f = Z[Ko + Ki-Z + K2(K; + Z)'] — 3, Zz xtaiv2 


f+iv2 
$sf<i 
(21) Ko = 10° X .18664 67623 69 K: = 10° X .13545 03944 219 
K, = .19531 25 = 0011001). K; = 10° X .13293 49397 97 


A (absolute error) + 2 x 10° 


Several of the above approximations listed below have been incorporated in an 
Elementary Function Subroutine Package for the IBM computers 704 and 709. 


Share Distribution 


Equation Machine Number Name 
(12), (14) 704, 709 510, 507 IB SIN1 
(16) 704, 709 510, 507 IB TANI 
(17) 704, 709 510, 507 IB TAN2 
(19), (20) 704, 709 571, 590 IB SIN2 
(21) 709 665 IB LOG3 


In conclusion, the author wishes to express his indebtedness to Dr. E. G. Kog- 
betliantz for his advice and guidance, and to Mr. F. S. Beckman, IBM, for the 
support of this project. 


IBM Corporation, New York, and 
College of the City of New York 
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418 N. REAM 


Note on “Approximation of Curves 
by Line Segments” 


By N. Ream 


The problem of obtaining a best fit of broken line segments to a curve over a 
given range has recently been investigated by Stone [1] who has prepared a general 
computer program to solve the least-squares equations. 

The problem arose previously in designing diode function-generators for analog 
computers [2], [3], [4]. If f(a) is the given curve and (wz , uy) is the range to be fitted 
by N segments, and if f(x) may be approximated by a parabola in each segment, 
then it may be shown [4] that the unweighted least-squares fit yields the following 
equation for the breakpoints uw; , --- , Unw-1: 


(1) [ 


0 


i {f"(a)}°"* dx = 2 fi tf" (x) }°* dz, 


and that the ordinate v; of each breakpoint is given by 


ee ie ere ft} ” (4, \ 0-2 on ‘en 0.4 F 
(2) 0) — flu) = ~ aah tru)" | "yr eny ae |. 
For f(z) = e ™ fitted over (0, 3), equations (1) and (2) become 


(3) 1 Ll e wens all (1 iu " as > 


i 
N 





(4) vj; — ge = a ars Dv gery’. 


Table 1 gives values of wu, and maximum error E,,,x computed from (3) and (4) for 
N = 2; Stone’s values are shown in parentheses. Z,,.x occurs atx = 0. The table 
also gives values of the r.m.s. error R which the least-squares analysis aims to 
minimize; R is computed from the formula 


(5) (uy — w)R® = (1/720N*) Lf 


uo 


{f” (a) }°" az | , 
which for the chosen function becomes 
(6) R = (6c) °° Eynax « 


The derivation of equations (1), (2), and (5) involves expanding f(x) in a Taylor 
series about the center of each segment and retaining the first three terms. Hence 
(i) the formulas are exact for a parabola—it follows immediately that the best fit 
to a parabola has equally-spaced breakpoints; (ii) the method fails where f”(x) = 0. 

It may be mentioned that if the “best fit” is required to minimize the maximum 
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TABLE 1 
f(x) = = fitted with 2 segments over (0,3) 
c wa Emax | R E 
0.1 1.454 | (1.385) 0.00166 | (0.0016) | 0.00215 | 0.00121 + 0 
0.2 1.410 | (1.400) 0.00593 | (0.0059) | 0.00541 | 0.00420 + 0 
0.3 1.366 | (1.360) 0.0119 (0.0119) | 0.00887 | 0.00821 + 1 
0.4 1.322 (1.316) 0.0189 (0.0189) | 0.0122 | 0.0127 +0 
0.5 1.278 | (1.276) 0.0265 (0.0265) | 0.0153 | 0.0174 +1 
0.6 1.236 (1.235) 0.0343 (0.0344) | 0.0181 | 0.0221 +1 
0.7 1.194 | (1.196) 0.0420 (0.0423) | 0.0205 | 0.0264 +2 
0.8 1.153 | (1.155) 0.0496 (0.0500) | 0.0226 | 0.0305 +3 
0.9 1.113 (1.116) 0.0568 (0.0574) | 0.0244 | 0.03438 +5 
1.0 1.074 (1.080) 0.0636 (0.0645) 0.0260 | 9.0377 +7 
1.2 1.001 (1.008) | 0.0758 (0.0774) | 0.0283 | 0.0435 + 13 
1.5 0.900 | (0.912) 0.0907 (0.0936) 0.0302 | 9.0500 + 26 





| 


error, the breakpoints are given by equations (1) with the index 0.4 replaced by 
0.5 and f”(z) replaced by its absolute value. The maximum error E£ is then given by 


(7) E= EB Es | f(z) | az 


For the function under discussion (7) becomes 


ithe. _ go i.bey2 

(8) E= aN? (1 é 

and the error 5£ in E due to the approximations used in deriving (8) may be shown 
to be given by 


(9) bE => hE’. 
Values of EZ and 6£ are included in the table. 


Battersea College of Technology 
London, 8. W. 11 
England 


1. H. Strong, ‘‘Approximation of curves by line segments,’’ Math. Comp., v. 15, 1961, p. 
40-47. 

2. H. Hamer, “Optimum linear-segment function generation,’’ Trans. Amer. Inst. Elec. 
Engrs., v. 75, 1956, p. 518-520. 

3. M. E. Fisner, ‘‘The optimum design of quarter-squares multipliers with segmented 
characteristics,”’ J. Sci. Instrum., v. 34, 1957, p. 312-316. 

4. N. Ream, “Approximation errors in diode function-generators,’’ J. Electronics Control, 
v. 7, 1959, p. 83-96. 





420 G. A. PAXSON 


The Compositeness Of The Thirteenth Fermat 
Number 


By G. A. Paxson 


Fermat numbers are numbers of the form F, = 2” + 1. As is well known, these 
numbers are prime for n = 0, 1, 2, 3, and 4. In fact, Fermat conjectured that they 
all were prime. Since then, however, factors of many F,, for n > 4 have been dis- 
covered [4]. In particular, factors are known for n = 5, 6, 9, 10, 11, 12, 15, 16, 18, 
23, 36, 38, 39, 55, 58, 63, 73, 77, 81, 117, 125, 144, 150, 207, 226, 228, 250, 267, 
268, 284, 316, 452, and 1945. On the other hand, Morehead and Western [1], [2] 
showed over fifty years ago that F; and F,; are composite. Thus, F,; was the first 
number whose character was unknown, although [4] it has no factors less than 2”. 

Since F,, is prime if and only if 


gira DR = —] (mod F,), 


a program was written for the IBM 7090 to check this criterion and determine the 
character of F;; . At the beginning of December 1960, the computation was made 
and the result shows that F;; is composite. 

The computation required 6 hours 17 minutes of machine time, including time 
to punch a number of intermediate residues. The result was checked by a rerun 
on a different day and all residues were identical. Moreover, as a check on the 
program, Morehead’s result [1] for F; , Morehead and Western’s result [2] for Fs , 
and Robinson’s result [3] for Fi were obtained. 

The calculation to determine the character of Fy, requires about 493 hours. This 
and the recalculation to check the result are each about half finished. However, 
lack of available machine time has temporarily suspended the computation. 

When machine time does become available, in addition to finishing the com- 
putation for Fy, , a search will be made to insure that Robinson’s list [4] of factors 
p of Fermat numbers contains all p < 2”. At present it contains all p < 2” and 
all p < 2” for which p = 1 (mod 2”). It should be noted that this last search 
could yield new factors only for 7 < n < 12 since all factors of F, are of the form 
k-2" 4 1. 

California Research Corporation 
Richmond, California 

1. J. C. Moreneap, ‘Note on Fermat’s numbers,” Bull. Amer. Math. Soc., v. 11, 1905, 
p. 543-545. 

YD. 2. J. C. Moreneap & A. E. Western, ‘‘Note on Fermat’s numbers,” ibid, v. 16, 1909, p. 
3. R.M. Rosrnson, ‘““Mersenne and Fermat numbers,’”’ Proc. Amer. Math. Soc., v. 5 1954, 
842-846. 


"4. R.M. Rosrnson, ‘‘A report on primes of the form k-2" + 1 and on factors of Fermat 
numbers,” ibid., v. 9, 1958, p. 673-681. 
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88[A-E, G, I, K, L, M, X].—Grantno A. Korn & Toeresa M. Korn, Mathematical 
Handbook for Scientists and Engineers, McGraw-Hill Book Co., New York, 
1961, xiv + 943 p., 24 cm. Price $20.00. 


Workers in every walk of mathematical life will gratefully welcome this latest 
addition to an illustrious series of Handbooks. In its variety and scope it may well 
be the largest collection of widely useful mathematical facts and data ever com- 
piled. It admirably fills a too-long existing gap among the handbooks available to 
workers in technical fields. As a “tool of the trade” its price is certainly reasonable, 
probably offering the “lowest cost per (mathematical) fact.” 

“This handbook is intended, first, as a comprehensive reference collection of 
mathematical definitions, theorems, and formulas for scientists, engineers, and 
students. Subjects of both undergraduate and graduate level are included. The 
omission of all proofs and the concise tabular presentation of related formulas have 
made it possible to incorporate a relatively large amount of reference material in 
one volume. 

“The handbook is, however, not intended for reference purposes alone; it at- 
tempts to present a connected survey of mathematical methods useful beyond 
specialized applications. Each chapter is arranged so as to permit rapid review of 
an entire mathematical subject,” and chapter introductions, notes, and cross- 
references interrelate the many topics “for a broad view of the entire field of mathe- 
matics.” 

To meet the requirements of different readers the material has been arranged 
at three levels: 

“1. The most important formulas and definitions have been collected in tables 
and boxed groups permitting rapid reference and review. 

“2. The main text presents, in large print, a concise, connected review of each 
subject. 

“3. More detailed discussions and advanced topics are presented in small 
print.” 

The following summary of the book’s twenty-one chapters and appendices gives 
a brief indication of the scope of the material. 

Chapters 1 through 5 review the basic college material on algebra, analytic 
geometry (plane and solid), elementary and advanced calculus, including Lebesgue 
and Stieltjes integrals, and vector analysis. Chapters 6, 7, and 8 cover curvilinear 
coordinates, functions of a complex variable, and Laplace and other integral trans- 
formations, respectively. Chapters 9 through 11 cover ordinary and partial differ- 
ential equations (including transform methods, method of characteristics), and 
maxima and minima, including the calculus of variations. 

Chapters 12 through 14 deal with various aspects of mathematical models. 
Chapter 12 introduces the elements of modern abstract language and covers con- 
cepts such as groups, fields, topological spaces, and Boolean algebras. Chapter 13 
deals with matrices, and quadratic and hermitian forms. Chapter 14 treats linear 
vector spaces and transformatiors, including matrix representation, eigenvalues, 
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and group representations. Chapter 15 handles the subject of linear integral equa- 
tions, boundary-value problems, and eigenvalue problems. Chapters 16 and 17 
give a good outline of the related subjects of tensor analysis and differential geom- 
etry. 

Chapters 18 and 19 recognize the increasing importance of statistical methods 
in many fields and provide over 100 pages devoted to probability and random 
processes, and mathematical statistics. The material is given in appealing detail, 
and the modern worker has the comfort of finding succinctly in a single source many 
of the not always easy-to-find formulas and results on such topics as multi-dimen- 
sional distributions, limit theorems, generalized Fourier analysis including correla- 
tion and power spectra, sampling distributions, and statistical estimation and 
testing of hypotheses. 

Chapter 20, on numerical calculations and finite differences, reviews the stand- 
ard methods and has a section on difference equations. Included are numerical 
methods for matrix inversion, eigenvalues, interpolation and approximation, 
ordinary and partial differential equations, and numerical harmonic analysis, 
among others. Chapter 21 is essentially a brief collection of formulas on the proper- 
ties of elementary and higher transcendental functions. 

A significant portion (about one-sixth) of the volume consists of six appendices 
as follows: formulas for plane figures and solids; plane and spherical trigonometry; 
permutations, combinations, and related topics; tables of Fourier expansions and 
Laplace-transform pairs; tables of indefinite and definite integrals; and twenty 
numerical tables. 

The book is rounded out with a glossary of symbols and notation showing where 
each item is explained, and a comprehensive index of almost thirty pages that en- 
ables the book to be used as a mathematical dictionary. 

The painstaking care with which each subject is organized is shown by effective 
use of summary tables and boxes. The box is a device in common use abroad, which 
might well be used more widely here. A small sampling of the tables and boxes will 
be helpful to the prospective user and give an insight into the valuable nature of 
the material: a table of formulas dealing with tangents, normals, and polars for 
each of the four classes of conic sections; a box showing various forms for the 
equation of a plane, and line, in both cartesian and vector notation; a table of 
properties of Fourier transforms—linearity, change of scale, shift, convolution, 
modulation, differentiation, Parseval’s theorem; tables of operations on scalar and 
vector point functions; tables relating to a wide variety of transformations and 
other properties of the various orthogonal curvilinear coordinate systems; a table 
of real and imaginary parts, zeros, and singularities of common functions; a graphi- 
cal set of sixty conformal mappings of regions in the complex plane; a table of 
definitions for different types of tensors; boxes with definitions of Riemann space 
and associated quantities such as covariant derivatives, Christoffel three-index 
symbols, and the first and second fundamental quadratic differential forms of a 
surface; a table explaining fourteen numerical parameters describing properties of 
one-dimensional probability distributions; boxed formulas for moments, character- 
istic and other generating functions for one-dimensional distributions, and for two- 
and more-dimensional probability and marginal distributions; tables for the many 
formulas and properties of a dozen or more discrete and continuous distributions 
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of most importance in applications, including the features “typical interpretation” 
and “approximations,” which are rarely presented in standard treatments in such 
useful form; tables of formulas for tests of hypotheses and confidence intervals for 
normal populations; a box for the unit-step functions, with accompanying sketches, 
and of relations involving the delta function and its “derivatives’’; and a table and 
sketches of various types of pulses and waveforms and their characteristics. 

In addition to the appended numerical tables there are several short tables in 
the chapter on numerical calculations: 5- to 7-place tables for Lagrange, Newton, 
Stirling, Bessel, Everett, and Steffensen interpolation, and tables for abscissas and 
weights for Gauss and Chebyshev quadrature formulas. 

While not detracting materially from the excellence of the book, mention of a 
few necessary corrections that were noted may be of help to the user. The pagina- 
tion for Chapter 6 should be corrected in the Table of Contents as follows: See- 
tions 6.4, 6.5, 6.6 begin on pages 170, 173, and 173, respectively; on page 112 the 
symbol ““— 0” is omitted under “lim” on the right-hand side of equation (4.6-43) ; 
on page 443, in the displayed equation at the top of the page, ““= 0” should be re- 
placed by “‘#0’’; on page 489, the second line under the last box, Sec. 5.10-3 should 
read Sec. 16.10-3; on page 566, the typography is confusing in the last line of the 
table owing to wrong size of type—the mathematical expressions should read 


1 a—1 —z/B 1 7 
BT(a) 2” e and Ta) TL (a); 
in the same table, the last two entries in the column “Characteristic function” are 
known in explicit form and should be given, namely, F(a; a + £8; it) (confluent 
hypergeometric function) and (1 — it)“; on page 570, in equation (18.8-29) the 
minus sign is omitted from the exponential, and in equation (18.8-30) the multi- 
plier (1/xa) is omitted from the expression for ¢.(x); on page 626 the reference in 
the heading of Sec. 19.8-2 should be “Sec. 18.12-2” instead of “See. 18.11-2”; on 
page 935, the index entry “Probability distribution” might well have included a 
reference to Sec. 18.8, as it has eight tables showing valuable information about the 
most important special distributions in statistics. 

It may be appropriate to mention several additional matters that may be of 
value in connection with any later edition. The subject of random numbers ap- 
parently is not included anywhere in the book, and it would seem that at least one 
of the most important modern works on this topic warrants mention either in 
Chapter 18, on Probability, or in Chapter 19, on Mathematical Statistics, namely, 
The Rand Corporation’s A Million Random Digits with 100,000 Normal Deviates, 
The Free Press Publications, Glencoe, Illinois, 1955. Also, some of the older refer- 
ences listed at the end of Chapter 19 should be replaced by their more modern 
versions; for example, Arkin and Colton is in a fourth, revised edition (1955), and 
P. G. Hoel is in a second edition (1954). In addition, the following works might be 
included as being very useful for reference and application purposes: 

Oscar Krisen Buros, Statistical Methodology Reviews, 1941-50, John Wiley & 
Sons, Inc., New York, 1951; 

M. G. Kendall & W. R. Buckland, A Dictionary of Statistical Terms, Hafner 
Publishing Co., New York, 1957; 

E. P. Adams & R. L. Hippisley, Smithsonian Mathematical Formulae and 
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Tables of Elliptic Functions, Smithsonian Miscellaneous Collections, Vol. 74, 
No. 1, Washington, D. C., 1922 (or later edition); 

Statistics Manuai, NAVORD Report 3369, Naval Ordnance Test Station, 
China Lake, California, 1955. 

Several minor points may be noted with regard to the numerical tables in Ap- 
pendix F’. To the eleven numerical constants listed should be added Euler’s constant, 
which occurs in a number of places in the text. There is space for increasing the 
number of decimal places shown to at least 10; this should be done to increase their 
usefulness. The typographical layout for several of the tables is hard on the eyes 
because little or no space is allowed between entries in adjacent columns. The 
columnar lines alone do not provide effective separation, so that the entries running 
across the page merge into one another. This applies to all or part of the tables for 
squares, integral sine and cosine, x? distribution, and F distribution. This can be 
remedied either by use of smaller type or by printing the tables along the length 
rather than the width of the page, as is done with some of the other tables, resulting 
in much greater legibility. 

Much of the material of the book is necessarily gathered from other sources. In 
a number of places, especially the figures, the source is cited from among the refer- 
ences at the end of the chapter. It would be helpful if such citation (admittedly 
laborious) could be done more systematically, as this could save a great deal of 
time and effort spent in searching through the listed references in order to follow 
up a particular theorem or development. 

As regards the physical aspects, one would wish that a book of such utility 
could be constructed in such a manner as to better be able to withstand the great 
amount of handling it is bound to receive, perhaps by being issued in the almost 
indestructible form achieved by the binders used in the tax and accounting services. 

Even with the minor shortcomings indicated here, this mathematical handbook 
is of such unique value that it can be unhesitatingly recommended for the intimate 
possession of everyone with a serious interest in the theory or application of virtually 
any aspect of mathematics. 


JuLius LIEBLEIN 
Applied Mathematics Laboratory 


David Taylor Model Basin 
Washington 7, D. C. 


89 [E, L].—L. N. Nosova, Tabliisy funkisit Tomsona i ikh pervykh proizvodnykh 
(Tables of Thomson Functions and their First Derivatives), Izdatel’stvo Akademii 
Nauk SSSR, Moscow, 1960, 422 p., 27 cm. Price 49 Rubles. 


This new addition to the series of tables prepared at the Computation Center 
of the Academy of Sciences, USSR, consists of two main tables. The first of these 
presents values of the Thomson (or Kelvin) functions ber x, bei x, ker x, kei x, and 
of their first derivatives to 7S for x = 0(.01)10. The second principal table gives 
values of modified functions consisting of the Kelvin functions of the first kind 
(ber, bei) and their first derivatives, each multiplied by e-*/V?, and the functions 
of the second kind (ker, kei) and their first derivatives, each multiplied by e*/V?. 
These data are also given to 78, for x = 10(.01)100. Corresponding values of e*/V? 
are tabulated in an adjoining column to 78. The entries throughout appear as 
7-digit integers multiplied by an appropriate power of 10. 





—_ ——— 











ca 


Se SV eee OP 


we 











REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 425 


The book closes with two smaller tables: the first consists of 78 values 
of exp (m-10°°/+/2) for m = 1(1)1000; the second gives exact values of t(1 — t)/2 
for t = 0(.001).500, for use in interpolation with second differences, which are given 
throughout the main tables. 

We are informed in the Preface that the underlying computations were per- 
formed on the electronic computer STRELA, using the well-known power series 
and asymptotic series for these functions. These details, as well as a discussion of the 
arrangement of the tables and their use, are presented in the Introduction. It is 
there stated that the tabular entries are each correct to within 0.6 of a unit in the 
last place shown. 

The reviewer carried out a partial check of the correctness of this statement 
concerning the accuracy of these tables, by comparing several entries in them with 
corresponding data in the recent tables of Lowell [1] which give values of the Kelvin 
functions and their first derivatives for z = 0(.01) 107.50 to between 9 and 14 sig- 
nificant figures. Comparison of the two tables revealed that the values given by 
Nosova for functions of the first kind and their derivatives are correct to within a 
unit in the last place for the range z = 0(.01)10. However, in the vicinity of zeros 
of the functions of the second kind and of their derivatives the Russian tabular 
values err by as much as 6 to 8 units, as, for example, kei zr when x = 8.24(.01)8.47 
and kei’(z) when z = 9.38(.01)9.43. The reviewer has verified, moreover, that 
ker’z is in error by 9 units when z = 7.16 and 7.18, and that ker’ 7.17 is too low 
by 86 units, which is explainable by virtue of the fact that this last tabular entry 
is only one-tenth its neighbors, and all three are subject to a nearly constant abso- 
lute error. 

This same difficulty in attaining the stated accuracy occurs in the second table 
in the book under review. Egregious examples of just a few of the large relative 
errors that were discovered occur in e*/V? kei x when x = 97.19 (tabular value too 
low by 335 units), in e-*/V? beiz when z = 98.30 (too high by 1027 units), and 
in e7/V2 ker x when xz = 99.41 (too high by 363 units). It should be stated here 
that the table of e*/V? was also checked at several places, and no errors were found. 

The arrangement of the material is very convenient, all functions of a given 
argument being found on facing pages. It is indeed unfortunate that the attrac- 
tiveness and convenience of these tables could not have been matched by acceptable 
accuracy. This accuracy could have been attained throughout by use of computer 
routines employing double-precision arithmetic, such as were used by Lowell. 

J. W. W. 

1. Herman H. Lowe 1, Tables of the Bessel-Kelvin Functions Ber, Bei, Ker, Kei, and their 
Derivatives for the Argument Range 0(.01)107.60, Technical Report R-32, National Aeronautics 
ss Administration, Washington, D. C., 1959. See Math. Comp. v. 14, 1960, p. 81 (Re- 
9o[H, S, X].—A. L. Logs, J. Ta. G. Overserx & P. H. Wrersema, The Electrical 

Double Layer Around a Spherical Colloid Particle, The Technology Press of 

M.1.T., Cambridge, Mass., 1961, 375 p., 26 em. Price $10.00. 


These tables give the numerical solution of the Poisson-Boltzmann equation 


--d(2dy\_ _ 4ne _ 2 ep\ _ ze 
r dr (+ #) iF In. 2+ exp ( se) nr 2 exp kT ] 


where y is the electric potential at radius r from the center of a charged, spherical 
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colloidal particle in an electrolyte; the local charge distributions and the free energy 
are also given. The electrolyte is characterized by the following parameters: «, the 
dielectric constant; z,, the valence of the positive ions; z_, the valence of the 
negative ions; n, , the concentration of positive ions far from the particle; n_ , the 
concentration of negative ions far from the particle; 7’, the absolute temperature. 

For the numerical computations reduced variables are introduced. In these new 
variables the Poisson-Boltzmann equation becomes 


dy _ exp(z_y) — exp(—zyy) 
dx? 2z_ x* 





and the boundary conditions are y = 0 at x = 0 andy = y at x = 2», where y 
is the reduced potential and x is the new independent variable. The local charge 
distributions and the free energy are represented by J,(z), J_(x), and F(x) where: 








ial “ 1- as r 
I,(z) =727 [ a dr; 

— a ’ e* we 1 . 
I(z) =z I "eas dr; 

eee #) z4(e°” — 1) — z(1 — €*+) 
(a) = # [5 (Gt) + ee ir 


The quantities x, y(z), J,(x), J_(x), and F(z) are tabulated for a variety of values 
of z,, 2 (1,1; 2,1; 3,1; 1,2; 1,3) and of 1/2» (from 0.1 to 20 in varying steps) 
and of yo (from 0.5 to 16 in varying steps). The values of y, J, , J_ and F are said 
to be accurate to four significant figures, except for a few cases where there is an 
error in the third figure. 

The tables include a forty-page discussion of the equation to be solved, the 
numerical methods and the results. 

These computations are said to be more extensive and more accurate than 
similar computations performed by N. E. Hoskin, Trans. Faraday Soc., v. 49, 
1953, p. 147. 

Luoyp D. Fospick 


Digital Computer Laboratory 
University of Illinois 
Urbana, Illinois 


91(H, X].—Loruar Cotiatz, The Numerical Treatment of Differential Equations, 
Third Edition, Translated by P. G. Williams from a supplemented version of 
the second German edition, Springer-Verlag, Berlin, 1960, xv + 568 p., 24 em. 
Price DM 98. 


The translation of Professor Collatz’s book into English will be welcomed by all 
those people who are in any way concerned with the numerical solution of differen- 
tial equations. No other single book on this subject contains such a vast amount of 
material. The following list of the chapter headings gives some idea of the range of 
topics covered. 

I. Mathematical Preliminaries and Some General Principles 
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II. Initial-Value Problems in Ordinary Differential Equations 

III. Boundary-Value Problems in Ordinary Differential Equations 

IV. Initial and Initial-Boundary-Value Problems in Partial Differential Equations 
V. Boundary-Value Problems in Partial Differential Equations 

VI. Integral and Functional Equations 

The appendix contains a number of tables giving the various numerical methods 
in handy tabular form for both ordinary and partial differential equations. 

This book is a translation of the second German edition with some differences. 
As the author states, “It differs in detail from the second edition in that through- 
out the book a large number of minor improvements, alterations and additions 
have been made, and numerous further references to the literature included; also 
new worked examples have been incorporated.” 

The book is large but by no means covers the subject completely, as the author 
is careful to point out in the preface. Professor Collatz also disclaims any attempt to 
make general critical comparisons of the various methods presented. This is to be 
regretted, since it decreases the value of the book to those persons who would be 
most likely to refer to it, namely the neophytes in this numerical field. Along this 
same line of criticism there is no mention made of the use of computers, either 
analog or digital, for the numerical solution of differential equations. This gives a 
new book a distinctly old-fashioned flavor. A specific example may be cited to 
illustrate the criticism. In Chapter II the author states that the Runge-Kutta 
and Adams methods are stable with respect to small random errors. However, he 
does not warn the reader as to which of the well-known methods are unstable. 
Thus, the uninitiated might be tempted to code an unstable method as a subroutine 
for a computer. 

On the positive side, the book can be recommended for its vast coverage, its 
many worked examples, and its close attention to error estimates. The translation 
is smooth and the printing is excellent. 

RicnHarp C. RoBEerts 


U. S. Naval Ordnance Laboratory 
White Oak, Silver Spring, Maryland 


92(H, X).—L. Derwipvus, Introduction a I’ Algébre Supérieure et au Calcul Nu- 
mérique Algébrique, Masson et Cie, Paris, 1957, 432 p., 25 em. Price 6000 fr. 
This author presents an interesting combination of pure theory and computa- 

tional methods for linear and non-linear algebra—that is, linear equations, eigen- 

values, roots of algebraic equations, etc. He concludes with an introduction to 
abstract algebra. 

The numerical methods are developed for use with desk computers and are 
amply illustrated throughout the text. A listing of the chapter headings will indi- 
cate the scope of the work. 

I. Mécanisation du calcul algébrique. Nombres complexes. 

II. Les déterminants et les systémes d’équations linéaires. 

III. Théorie générale des polynémes et des fractions d’une indéterminée. 
IV. Elimination et systémes d’équations algébriques. 
V. Résolution numérique des équations. 
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VI. Substitutions linéaires, formes quadratiques et transformations rationnelles. 
VII. Calcul matriciel 


VIII. Equations dont les racines sont dans un cercle ou un demi-plan. Critéres de 
stabilité. 
IX. Notion sur les groupes et sur l’algébre abstraite. 
Appendice Sur les determinants de Hurwitz et la séparation des racines complexes 
des équations & coefficients réels. 


E. I. 


93[H, X)—Minorvu Urase, Hrroxt YANAGIwARA & YOSHITANE SHINOHARA, 
“Periodic solutions of van der Pol’s equation with damping coefficient } = 
2 ~ 10,” reprinted from the J. Sci. Hiroshima Univ., Ser. A, v. 23, No. 3, 
March 1960. 


The periodic solution of van der Pol’s equation 


dy” _ 2) dx r 
de ACL xJar+x=6 


is tabulated for A = 2, 3, 4, 5, 6, 8, 10. For each \ a four-decimal-place listing of 
the function x(t) and the function 


dx 
= < 
at for A~<4 


y(t) = 
% | for A>5 
dt thas 


is given for the range 7,(a) S ¢ S T2(a), where a is the initial positive amplitude 
of the periodic solution normalized so that at t = 0, x = 0; and where 7;(a) is 
the smallest positive time at which x = 0, while 7;(a) is the largest negative time 
at which x = 0. 

Since the periodic solution corresponds to a closed curve in the (x, y) plane 
which is symmetric with respect to the origin, the above tabulation is sufficient. 


For \ > 5, an additional three-decimal-place tabulation of x is given. 
The interval size in ¢ depends on \ and on the value of ¢ as in the following table: 














r t>0o0 t<0 
2 .05 .025 
3,°°°,8 .025 .0125 | 
0.0(.0125)0.2 
10 0.2(.025)9.0 .00625 
9.0(.0125)9.25 








Hach table includes a listing of the same quantities at = 7,(a) andt = T.(a). 
Furthermore, four-decimal values of the amplitude a, the period » = 
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. 


2(T:(a) — T,(a)], and the characteristic exponent h are given for each X. If we set 


h(t) =X [ (1 — x’) dt, then h = h(w)/o. 


x) and the curve x(t) (including 
\ = 0, 1). An additional graph depicts a, w, and h as functions of ) in the interval 
[O, 10}. 


For each ) there is a plot of the hodograph (x, 


E. I. 


94[M, P, S].—R. L. Murray & L. A. Minx, Tables of Series Coefficients for Burnup 
Functions, Bulletin No. 71, Department of Engineering Research, N. C. State 
College, Raleigh, N. C., May 1959, 82 p., 28 cm. Price $1.50. 


In a certain model, calculation of nuclear reactor properties under long-term 
operation requires the evaluation of 


ne 1 6.5/2 ; I 2 5-50 | 
o=@ F ~5 /, (cos x)‘ dx iam | z[Jo(x)}° dz |, 


_ gD b igh 1 fo . 
at ape Ag 
¢g ¢g 

















and some other combinations of yg. Here jo is the smallest positive zero of Jo(x). 
All functions are tabulated for the range 1 = 0(1)4, 4; , dg = 0.50(0.05)0.60(0.02)1.0. 
Ao and dp are given to 7D; A; to 6D, and the remaining functions not listed here 
are given with less accuracy. Similar tables are given for the function 


Ao = aD ie x (= *) dx 
© Gry? Jo x 
The method of computation is not explained, nor is jo defined. We infer the defini- 
tion of jo from physical considerations corroborated by numerical evaluation. Spot 
checks indicate the entries are accurate to the number of places given. 
Y. L. 





9o5[W, X].—Russextt L. Acxorr, Editor, Progress in Operations Research, Vol. 1, 
John Wiley & Sons, Inc., New York, 1961, 505 p., 23 em. Price $11.50. 


Each chapter of this book is written by different authors. It treats recent prog- 
ress in some of the methodological fields of operations research, such as linear pro- 
gramming, in an outstanding manner, scantily discussing progress in others, such 
as queuing theory. 

The introductory chapter, written by the editor of the book, is excellent. He 
points out that in other well-established scientific fields one is not as concerned 
about definitions as those in operations research have been, that this field is now 
accepted and has acquired the confidence of workers in other fields, and that, as a 
result, there is less craving for definitions. 

An interesting chapter by Churchman on contributions to decision and value 
theory then follows. Hanssmann’s chapter on inventory theory leaves much to be 
desired and is not saved even by attempting to justify the presentation in an opera- 
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tions research framework. The use of parallel and series stations leaves many prob- 
lems untouched. The stochastic nature of the field does not come through satisfac- 
torily. A chapter on mathematical programming follows, in which Arnoff and 
Sengupta give a superb account of progress in programming except for non-linear 
programming, in which there have been several recent contributions, such as that 
of Zoutendijk [1]. A remark at the top of page 176 regarding the unavailability of 
work on sensitivity is inaccurate. This reviewer has proved in the 1959 paper re- 
ferred to on page 209 that at a solution vertex the objective function (in the cus- 
tomary notation) has the following sensitivity to a;; 


av _ _ .%0 _ oVav 
1” oh dc; 


04;; 
where z;’ and y;’ are the solutions to the primal and the dual, respectively. A read- 
able and very useful account of dynamic programming, including adaptive proc- 
esses, is then given by Dreyfus. Chapter 6 by Morse deals with Markov and queuing 
processes. Sisson studies sequencing theory in chapter 7, and a variety of very 
useful replacement models, developed by a number of individuals, are treated by 
Dean in the next chapter. In another chapter, Morgenthaler describes simulation 
and Monte Carlo in a manner which provides useful guide-lines for application. 
Thomas, well known for his contributions to game theory, treats the subject in 
chapter 10 in an interesting style which utilizes historical ideas on the subject. The 
presence of these two chapters clarifies in the mind of the reader differences between 
simulation and gaming. Magee and Ernst examine the future of operations research 
in chapter 11 and point out the need of quantitative models of human behavior, 
marketing, interaction of men with men and with machines, organization, and in- 
formation; and they call for a better grasp of risk and competition. They point out 


that operations research is far from mature and has promise. This book is recom- 
mended reading. 





Tuomas L. Saaty 
Office of Naval Research 
Navy Department 
Washington, D. C. 


. G. ZouTenpisK, ‘Maximizing a function in a convex region,” J. Roy. Statist. Soc., 
1B. 1959, p. 338-355. 





Methods of Operations Research, McGraw- 
Hill Book Co., Inc., New York, 1959, xi + 421 p., 24 em. Price $10.00. 


This book is about some selected mathematical methods of operations research, 
but it offers both more and less than what its title may suggest to some readers. 
Though full of mathematical results, this book is not a cycle of “lemma, theorem, 
proof, and corollary.” Though it has many problems, and though Saaty is deeply 
concerned with principles of solution, this book is not a ‘problem manual.’”’ Despite 
the fact that there are many illustrative applications, this is far from being either 
a “cook book” or a collection of case studies. 

For the unifying thread of the book, one must look to the exuberant creativity 
of Saaty himself. He sets the tone of the book in the preface with the statement 
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that begins, ‘““‘We wish to warn the reader at the outset that even though we may 
attempt partly to inform him and partly to stimulate him .. .” It is the effort to 
stimulate the reader, to jolt him from his mental ruts, that is most striking. Saaty 
continually reminds the reader of the elemental role of creativity, that formal proof 
is not all of mathematics, nor is mathematics all of operations research. Operations 
research, incidentally, according to Saaty’s preferred definition, “is the art of giving 
bad answers to problems to which otherwise worse answers are given.” 

The preface states nominal reader requirements of “‘a course in calculus, with 
some elements of advanced calculus and rudimentary knowledge of matrix theory,” 
but it also admits to “compromise in the presentation” as a response to an antici- 
pated “‘variety of background among the readers.’”’ As the development unfolds, it 
becomes clear that Saaty hesitates neither to emphasize an elementary result that 
he finds suggestive nor to introduce a more advanced result that he finds particu- 
larly intriguing. 

The outline of the book is much what one might expect from the title. The in- 
troductory first chapter presents history and concepts of operations research. Then, 
Part 1 gives three chapters, also largely introductory in nature, that range rather 
widely over topics related to the scientific method as an approach to truth, to mathe- 
matics and logic as approaches to validity, and to some elementary classical methods 
useful in the formulation of mathematical models. As the first of the two major 
sections, Part 2 has three chapters on the subjects, respectively, of optimization, 
linear and quadratic programming, and the theory of games. In Part 3, the other 
major section, there are four chapters that dispose of basic probability, applications 
of probability, fundamental statistics, and queuing theory. The lone chapter of 
Part 4 concludes the book with an unusual essay on creativity. Each of the four 
parts gives a collection of interesting problems, and each chapter ends with a set of 
valuable references. 

The logic of the outline does not always withstand successfully Saaty’s attempts 
to stimulate the reader. Sometimes Saaty seems to give way to his own effervescence, 
getting ahead of his story by appealing to concepts or methods that are stated 
clearly only in a later chapter. The resulting unevenness of presentation would be a 
fault in a conventional text with more staid objectives. But conventional standards 
hardly apply to a book that draws upon such diverse topics as James Thurber’s 
rooster, the trial of Madeline Smith for the murder of her lover, ‘““The Critique of 
Pure Reason,” “The Prince” of Machiavelli, and ‘“The Three Princes of Serendip.” 
Even the typographical errors, which are of at least the usual frequency, may be 
viewed, generously, as useful stimuli to the reader’s alertness. 

There is only one reference to computers (“‘Computers are often used in gaming 
methods... They have the advantage of speed and economy and can be con- 
trolled and relied upon.”’). There are several passages that touch on topics of nu- 
merical analysis. What the book conveys in fullest measure, however, is the spirit 
of mathematics and a passion for the creative solution of problems, wherever they 
may arise. 

Ciayton THOMAS 


U. 8. Air Force 
Washington 25, D. C. 
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97(X, Z|.—_UNESCO, Information Processing, Proceedings of the International Con- 
ference on Information Processing, R. Oldenbourg, Miinchen & Butterworths, 
London, 1960. Distributed by International Documents Service, Columbia 
University Press, 520 p., 30 cm. Price $25.00. 


This volume publishes the proceedings of the first full-scale international con- 
ference on information processing by the use of modern digital high-speed calcu- 
lators. The conference was sponsored by Unesco, and was held at the Sorbonne and 
at Unesco House in Paris during 15-20 June, 1959. It was attended by nearly 2000 
participants from 37 countries. Fifty-nine papers were presented at eleven plenary 
sessions; in addition, approximately sixty short lectures were given at twelve sym- 
posia. The volume is divided in seven main chapters covering the following areas: 

Chapter I—Methods of Digital Computing 

Chapter II—Common Symbolic Language for Computers 

Chapter III.—Automatic Translation of Languages 

Chapter IV.—Pattern Recognition and Machine Learning 

Chapter V.—Logical Design of Computers 

Chapter VI.—Special Session on Computer Techniques of the Future 

Chapter VII.—Miscellaneous Topics 
In addition, there are published introductory or closing remarks by the Editor, 
S. de Picciotto; Réne Maheu, Director General of Unesco; Howard H. Aiken, Presi- 
dent of the Conference; André Danjon, President of the Association Frangaise 
de Calcul; Pierre Auger, Secretary General of the Conference; and Hughes Vinel, 
representative of the Department of State of France in the field of scientific research. 

The list of authors and participants includes some of the best-known names in 
the field of high-speed calculators and their application. Thus, although the papers 
are not uniformly excellent, there is no doubt that the material contained in the 
volume constitutes the most comprehensive compilation of knowledge in certain 
areas of information processing available at the time of the meeting. The talks are 
especially lucid, and the discussions extremely helpful in clarifying many points. 
Also, the summaries, presented by the chairmen of the various sessions, are well 
done and easily readable. Of special interest are the frank and clear exchanges of 
technical information between the U. 8. specialists in the field of MT (machine 
translation) and their counterparts in the USSR. D. Panov sets the tone for these 
meetings in his introductory statement, which is punctuated with humorous remarks 
(almost “wisecracking” in character) in the best American style. Thus he divides 
the work in the field of MT in four stages, (1) “Talking” about future accomplish- 
ments, (2) “Complacency,” when the algorithm has finally been constructed and 
we say how good it would be if only it worked, (3) ‘““Enthusiasm,”’ when the algo- 
rithm works but still has many difficulties, and (4) ‘The final stage—more talking.”’ 

Because of the large number of technical papers, and the wealth of material 
contained in the volume, it is not possible to give special attention to any number 
of selected papers in this review. The accession of this historic volume is a ‘“‘must”’ 
for every modern technical library. 

The sponsorship of this conference by the United Nations and the success which 
it enjoyed tend to emphasize the significance of this field of research in the modern 
world of great technological advances. There is no doubt that it will serve as the 
forerunner for many other such meetings to be held in the future. 


H. P. 
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98[Z].—Ivan Fores, Computer Logic: The Functional Design of Digital Computers, 

Prentice-Hall, Inc., New Jersey, 1960, xii + 458 p., 24 em. Price $12.00. 

The author has attempted to write an all embracing book on computers for the 
reader who has some scientific training, but little engineering or mathematical train- 
ing. In the reviewer’s opinion, the author “talks down”’ to his reader and interlaces 
his discourse with many irrelevant comments; e.g., “The concept of the individual 
is one of far reaching consequences. Would you expect this idea to have ramifica- 
tions in biology, law and ethics? . . .”’ (from section 7.1 of Chapter 7, Number Sys- 
tems and Counting). 

The author’s style is unfortunate in another respect: he often uses a term or an 
expression before he has defined it. In some cases the definition can only be found 
in the glossary at the end of the book. The reader is not helped by the fact that 
certain terms used in the text and partially defined therein are not listed in the 
index; e.g., the term “overflow’’. 

It is unfortunate that the book is not error free. Thus the unwary reader will 
have difficulty with the following passage from page 100 where the base 12 numbers 
are being discussed: ‘““To represent the quantity thirty-two, we would count out 
our bundle of a dozen twice and have seven (sic) units left over. Thus thirty-two 
will be represented . . . (by) . . . (27): .” 

The first half of the book is intended to give “‘a bird’s-eye view from the air . . . 
to see how the computer fits into the over-all system of scientific investigation and 
business enterprise.” The second half discusses the logical design of various basic 
units of a computer and the synthesis of larger functional units of such a machine. 

The sequence in which the author takes up various topics seems somewhat 
strange. Thus the chapter on Number Systems and Counting follows the one on 
Machine Arithmetic. In the latter chapter, addition is done by use of addition 
tables but the reasons for using these instead of counting techniques is never 
discussed. 

The list of chapter headings follows: 


Chapter One Introduction 

Chapter Two First Principles and Definitions 

Chapter Three Specifying the Computer For the Problem 
Chapter Four The Flow and Control of Information 
Chapter Five Coding 

Chapter Six Machine Arithmetic 


Chapter Seven Number Systems and Counting 
Chapter Eight Machine Languages 

Chapter Nine Logic 

Chapter Ten Logical Construction 

Chapter Eleven Functional Units 

Chapter Twelve The Logic of Arithmetic 

Chapter Thirteen Memory Devices and Their Logic 
Chapter Fourteen The Control Unit 

Chapter Fifteen Input and Output Equipment 
Chapter Sixteen A Problem 
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99[Z|.—F ritz Reutrer, Die nomographische Darstellung von Funktionen einer kom- 
plexen Verdnderlichen und damit in Zusammenhang stehende Fragen der prak- 
tischen Mathematik, DK 518.3:517.53. Forschungsberichte des Landes 
Nordrhein-Westfalen, herausgegeben durch das Kultusministerium, Nr. 912, 
Westdeutscher Verlag, K6ln and Opladen, 1960, 123 p., 30 em. Price DM 35.40. 


This volume is primarily a collection of material given by the author in a series 
of previously published papers. The first part of the book deals with methods to 
represent the analytic function w(z) = u(z,y) + iv(z,y) by a nomogram with two 
parallel nonlinear scales and two “sliding” curves. For example, a reading-line 
tangent to the curves u(z,y) = c and v(z,y) = c2, where c, and c are constants, 
intersects the parallel scaies in points from which zx and y may be read. Variations 
of this technique are also considered. For instance, the parallel scales may represent 
u and v while the sliding curves represent x and y. 

Examples include the elementary functions w = 2, z = w'*®, w = sinz and 
w = Inz. Considerable attention is devoted to the Jacobian and Weierstrassian 
elliptic functions. The 30 illustrative nomograms which make up the second part of 
the volume are very neatly drawn and easily read. The author also includes on each 
figure at least two numerical problems so that the reader can quickly become 
familiar with their use. 


Y. L. L. 


100(Z|.—Marsnat H. Wruset, A Primer of Programming for Digital Computers, 
McGraw-Hill Book Co., Inc., New York, 1959, xv + 230 p., 24 em. Price $7.50. 


As a good professor presents a new topic through ideas which are already fa- 
miliar to his students, Professor Wrubel, in the introductory chapter of A Primer 
of Programming for Digital Computers, smoothly leads his reader to a basic under- 
standing of the electronic digital computer and of the nature of programming. The 
potential programmer is cautioned not to attribute superhuman powers to the 
computer, and he is advised to consider carefully the value of a computer solution 
to each individual problem. 

Written for “scientists, engineers and their students who are planning to use 
computers as tools of research,” the primer is divided into two sections: Elementary 
Programming and Advanced Programming. The first section instructs the novice 
in the elements of programming as illustrated by the Bell Laboratories interpretive 
language. The second section trains the more experienced programmer in advanced 
techniques and in the use of a machine (IBM 650) language. Chapter subheadings 
in both sections follow the same sequence, enabling the reader to easily apply his 
basic knowledge to the more complex programming concepts. 

In his approach to the elements of programming Professor Wrubel defines and 
builds on each component from the basic digit through the word and the language, 
until he completes his construction of the computer program. Arithmetic, logical, 
and input-output control instructions are taught with clarity. Taking advantage of 
his readers’ familia. ity with mathematical notation, the author employs symbolism 
freely in explaining such aspects of programming as conditional transfers, address 
modification and looping, and subroutines. Professor Wrubel also describes and 
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stresses the importance of drawing the flow diagram, of code-checking, and of pre- 
paring the written report of each problem. 

The presentation of the Bell code and of the elementary programming concepts 
is lucid, and provides adequate information to enable the novice to program his 
problem. In this part of the text, particularly, most of the reader’s questions are 
anticipated, and it is obvious that the author speaks with a great deal of practical 
experience. The chapter of instructions for problem testing probably should be 
supplemented by IBM publications related to this phase of programming. 

In the final chapter of this section, Professor Wrubel speaks briefly on the sub- 
ject of automatic programming, and then immediately re-introduces the elementary 
concepts in the language of FORTRANSIT, an automatic programming system for 
the IBM 650. This chapter might have been more meaningful to the reader had 
the author’s commentary on the general nature of automatic programming been 
more fully developed. 

The Advanced Programming section of the text treats the basic functions of 
the computer through the machine’s own code of instructions and through SOAP, 
the Symbolic Optimal Assembly Program. The reader becomes acquainted with the 
use of error conditions for program control, with double-precision arithmetic, with 
the scaling of variables and constants, and with many more advanced topics. 

The primer includes many practice problems and a useful glossary of program- 
ming terms. 

Although the primer’s discussions are primarily concerned with programming 
for the IBM 650, they cover the concepts of programming in such a way as to be 
valuable to all newcomers to the field. 

Rrra Horsetr Burns 
Research Computing Center 


Indiana University 
Bloomington, Indiana 





NOTES 


SIAM—1U.S.S.R. Collaborative Publication 


On February 16, 1961, The Society for Industrial and Applied Mathematics 
(SIAM) announced publication of the English edition of Problems of Continuum Me- 
chanics, a volume dedicated to Academician N. I. Muskhelishvili on his seventieth 
birthday. On this same day, the Academy of Sciences of the U.S.8.R. announced 
publication of the Russian edition, and these companion volumes have been pre- 
sented to Muskhelishvili. This event is of interest, not only because the volume 
contains the recent work of leaders in the field of continuum mechanics, but also 
because it is the result of a collaborative publication effort between the Academy 
of Sciences of the U.S.S.R. and SIAM in the United States. 

SIAM obtained the cooperation of the Academy of Sciences of the U.S.S.R. to 
prepare an English edition which could be published simultaneously with the 
Russian edition. To enable this to be done, the Academy of Sciences supplied SIAM 
with manuscripts of contributed papers in sufficient time to permit simultaneous 
bilingual publication. 

The editorial board of the Russian edition is as follows: A. V. Bitsadve, G. Iu. 
Dzhanelidve, 8. A. Khristianovich, M. A. Lavrent’ev (Editor-in-Chief), A. I. 
Lur’e, G. F. Manvzhavidze, G. K. Mikhailov (Editor), L. I. Sedov, D. I. Sherman, 
S. L. Sobolev, V. V. Sokolovskii, I. N. Vekua (Assistant Editor). 

The editorial board for the English edition consists of I. E. Block and J. R. M. 
Radok. 


New Journal 


BIT, Nordisk Tidskrift for Informations Behandling, is a new quarterly journal 
in the computer field, supported by the Scandinavian countries, including Finland 
and Iceland, and published by the Societies for Information Processing in these five 
countries. 

BIT will contain scientific articles on numerical analysis, computer-technical 
problems, programming, operations research, and data-processing, and will publish 
informative articles on important new subjects. It is also intended to publish 
algorithms written in ALGOL. 

Scientific papers of the kind mentioned above are welcomed from foreign coun- 
tries. Further information concerning this new journal may be obtained from Dr. 
Carl-Erik Fréberg, Editor, Lunds Universitet, Solvegatan 14, Lund, Sweden. 


CORRIGENDA 


Harry WEINGARTEN & A. R. D1 Donato, “‘A table of generalized circular error,’’ Math. 
Comp., v. 15, 1961, p. 169-173. 

In the table of values of the generalized circular probable error, which appears on p. 170, 
the entry corresponding to c = .25, P = .55 should read 0.80030 instead of 0.80039. 

Hersert E, Sauzer & Norman Levine, “Table of a Weierstrass continuous non-differen- 
tiable function,’ Math. Comp. v. 15, 1961, p. 120-130. 

On the bottom line, p. 125, for pe 1/m3 ~ 1/2n? ~ O(h) read pr 1/m* ~ 1/2n? ~ oth). 
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VOLUME XV 
AUTHOR INDEX 
Papers and Technical Notes 


Short papers and notes are marked (N) in this index 


Author Title Page 
Arruurs, A. M. & Expansion of spherical Bessel functions in a series of 
McCarro., R. preg I Ro Simpy 159 
BarakatT, RICHARD Evaluation of the incomplete gamma function of imaginary 
argument by Chebyshev polynomials......... 7 
Best, Gitpert C. Two theorem tables of matrix algebra....................... 19 
Boersma, J. Two formulas relating to elliptic integrals of the third 
RE 5 eRe SS thre 44 . 296 
Burcuer, J. C. A partition test for pseudo-random numbers (N)............ 198 
Cooter, J. W. An improved eigenvalue corrector formula for solving the 
Schrédinger equation for central fields. . Serpette’, 
CorrineTton, Murian 8. Applications of the complex exponential integral. coef berries 1 
Dr Donato, A. R. & Integration of the geuveral bivariate Gaussian distribution over 
JARNAGIN, M. P. hl aes gpa gear 9 oe et I dete aye gdp 375 
Dranorr, J. 8. An eigenvalue problem arising in mass and heat transfer 
al ple pill cea aeamety » 5 ree myn emt «- 403 
Fre.ps, Jerry L. & Expansions of hypergeometric functions in hypergeometric 
Wimp, Jer ge tills a lie yg de paipegheny ong Megara be OPTI Ts 
Fruuer, L. & Stability analysis and integration of the viscous equations 
Lupuorr, H. F. thts ben drab ppabdend tire | ot ieee ee nyt yy srry ye 261 
FREUDENSTEIN, Bi-variate, rectangular, optimum-interval interpolation (N).. 288 
FERDINAND 
Gavutscut, WALTER Recursive computation of the repeated integrals of the error 
NS ds SP enseidutir a na apne aint alam eae . 227 


GREENBERGER, Martin’ An a prior: determination of serial correlation in ‘computer 
generated random numbers................... 


Hasetier, G. J. & Symmetric successive overrelaxation in solving diffusion dif- 
Wacuspress, E. L. EB eer Reb nt sas ed cenagesnes gues ved 356 
Hase.erove, C. B. A method for numerical integration....................... . 323 
IsemonceRr, K. R. Complete factorization of 2-1 (N).................... .. 295 
Jonss, J. G. On the numerical solution of convolution integral equations 
and systems of such equations.............. banecnd ae 
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